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ABSTRACT 


The report presents the analysis background and numerical details 
of a computer code for the transient thermal response and ablation of non- 
charring heat shields and nose-tips. The controlling sub-code is the 
in-depth transient heat conduction routine? it treats two-dimensional axi- 
symmetric bodies using the explicit finite-difference method. It is coupled 
to a completely general thermochemistry sub-code for the surface state and 
ablation computation. This routine allows for full equilibrium considering 
a large number of species, and can also account for any number of non- 
equilibrium effects. The code also features automatic coupled calculations 
of the pressure around the body, using a variation of modified Newtonian 
relations coupled to Prandtl-Meyer expansions. The subsonic region pressure 
description includes empirical corrections necessary on very blunt bodies. 

In addition, the code automatically computes the boundary layer edge static 
enthalpy, the recovery enthalpy, and the convective transfer coefficient, 
employing appropriate computed property values for any edge gas. 
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SECTION 1 


GENERAL INTRODUCTION 


1 . 1 PROBLEM STATEMENT 
1-1-1 Background 

The computation of the thermal and dynamic history of a vehicle entering 
a planetary atmosphere has a number of important features, which may be 
crudely summarized as follows: 

I. Trajectory analysis 
II- Inviscid flow field calculation 

III. Boundary layer calculation 

IV. Interactions between vehicle surface and boundary layer 

V. In-depth response 

Complete calculation of vehicle response over a period of time would involve 
a coupled calculation of all five phases, each computation phase providing, 
in effect, the time dependent boundary conditions for the adjacent phases. 

A recent survey (Ref. 1) has carefully examined the current computa- 
tional state-of-the-art for each phase and for the coupling between phases. 
Major deficiencies were identified in Phase III, boundary layer calculations, 
where general solution procedures for reacting multicomponent boundary layers 
did not exist, in Phase IV, where sufficiently general chemical routines did 
not exist, and in Phase V, where multidimensional in-depth solution programs 
did not exist. In addition, no III-IV-V coupled programs existed. Programs 
with IV- V coupling existed, but were limited to one space dimension for the 
in-depth response. 

Efforts to upgrade the computational state-of-the-art will probably 
not attempt to couple together large "complete" computer routines for all of 
the five phases? the resulting computer code would be unwieldly, uneconomical, 
and unreliable. Since, in general, the coupling aspects of the predictive 
problem are more interesting and important than the degree of computational 
detail in the individual phases, a more judicious computational approach seeks 
to couple detailed calculations of only one or perhaps two problem phases to 
less detailed approximate or correlation calculations of several of the other 
phases. The resulting computer code would include the important coupling 
effects without becoming excessively large, and would still offer good detail 
in at least one phase of the total problem. 



1.1.2 Present Computer Code 

The work reported here exemplifies this latter approach. The code 
developed under this program includes large subroutines for the detailed 
computation of the transient in-depth temperature response of an ablating 
heat shield (Phase V in the scheme above) and also of the thermochemical 
interactions between the environment gas and the heated surface of the heat 
shield material (Phase IV) . The program features more simplified approaches 
for the boundary layer transport calculations (Phase III) , for which a con- 
vective transport coefficient approach was adopted, and for the inviscid 
flow field calculations (Phase II) , for which empirical correlation and various 
modifications of Newtonian pressure distributions are used to find the pressure 
around the ablating body. These combine with the isentropic expansion assump- 
tion to define the boundary layer edge state. 

The completed code may thus be termed a "II -HI - IV - V coupled code", 
in which the routines for phases IV and V are relatively complete, general, 
and large, while those for phases II and III are relatively small and simpli- 
fied. The nature of the individual phase routines may be summarized as 
follows : 


Code Characteristics 

2-D Axisymmetric Finite Difference 
Transient Heat Conduction Code, 
Explicit, Multi-material Capability, 
General Geometry but some restric- 
tions on grid layout make use on 
very sharp bodies difficult, non- 
charring materials, some types of 
anisotropy considered 

Phase IV - Surface Response General thermochemical state routine 

allows for ablation of any one 
material in any environment, con- 
siders full equilibrium plus any 
number of heterogeneous and homo- 
geneous (but surface area scaled) 
kinetically controlled reactions; 
identifies thermo chemically control- 
ling surface material, detects failing 
or melting compounds 


Phase 

Phase V - In-Depth Response 
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Phase 


Code Characteristics 


Phase III - Boundary Layer 


Phase II - Flow Field 


Phase I - Trajectory 


Transfer Coefficient routine based 
on energy integral technique, 
laminar only at present time, 
boundary layer radiation not evaluated 

Modified Newtonian pressure for 
sharp bodies , blending into empirical 
pressure correlations for blunt 
bodies combined with Prandtl-Meyer 
expansions for supersonic regions; 
isentropic expansion used for boundary 
layer edge and recovery state cal- 
culations , flow field radiation 
not evaluated 

Not included in program, time 
histories of stagnation pressure 
and enthalpy to be input 


1.1.3 Purpose of This Report 

This report presents a description of the coupled code. It describes 
the general nature of the code and the flow of information from one computation 
phase to the other. The report also describes the physical models used in 
each computational phase, and outlines in some detail the numerical analysis 
used in each phase. The report does not include a users guide; this has been 
published separately (Reference 2) . 
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SECTION 2 


OVERVIEW OF COMPUTER PROGRAM AND REPORT 
Figure 1 illustrates the general physical problem considered by the 


code . 


Physics 

Main stream 


Boundary layer edge 



Computation 

Calculation of radiation 
flux to wall 


Determination of 
boundary layer edge 
state 


:count for boundary layer trans- 
port of mass and energy to and 
from surface 

Determine chemical state at sur- 
face and perform energy balance 
coupling to in-depth s lution 

In-depth thermal response 
calculation 


Figure 1. Sketch of Problem Considered 

In the coupled code, the routine for the in-depth thermal response 
calculations (Phase V) functions as the main, controlling routine. It is 
referred to as MACABRE, after "Material Ablation with Chemically Active 
Boundary Layers in Re-Entry". The coupling between Phases V and IV effected 
through the surface energy balance calculation, which is performed by a 
separate subroutine, SURFB. This energy calculation naturally requires (among 
other things) information about the enthalpy of the boundary layer gases 
adjacent to the wall; this information is called for by SURFB when needed 
and obtained from a package of complex thermochemistry subroutines controlled 
by subroutine EQUIL. In addition, SURFB requires values for a convective 
transfer coefficient, obtained as needed from subroutine RUCH , for the local 
pressure, obtained as needed from subroutine RUNLP , and for the local recovery 
enthalpy, obtained as needed from subroutine PARBL . Calls to PARBL , RUNLP, 
and RUCH constitute IV - III - II coupling, filling out the coupling require- 
ments of the code. 
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Sections 3 through 9 below describe the various computational phases 
and the coupling links between them. This necessarily involves much detail. 
Section 10 then presents a recapitulation of this introductory section and 
some overall flow charts of the computer code operations which should illus- 
trate more clearly all of the important features of the program. 



SECTION 3 


IN-DEPTH TRANSIENT THERMAL CALCULATION 

3.1 GENERAL REMARKS 

The in-depth transient temperature calculation serves as the control- 
ling logic unit for the entire coupled computer code. (It also consumes most 
of the computing time.) This calculation is of the familiar finite differ- 
ence, lumped-capacity type. The following sections describe the nodal grid 
employed for this calculation and the basic equations employed. 

3.2 NODAL LAYOUT AND GEOMETRY 
3.2.1 General Pattern 

The geometric shape considered, illustrated in Figure 2, is imagined 
to be divided up into a grid pattern in the customary manner for finite dif- 
ference computation schemes. The area within each grid "box" will be termed 
a "node"; this is in contrast to terminology which calls each corner of the 
grid network a node. The thermal capacity of each node is imagined to be 
lumped at a single point within the box; this point will be termed the nodal 
center. It will be convenient to assume for the moment that the nodal cen- 
ter may be located anywhere within the nodal box. Strategies for optimizing 
the location of the nodal centers will be discussed later. 

For convenience, the nodes are imagined to be quadrilaterals, so that 
the entire nodal network is an assemblage of quadrilaterals. For bookkeeping 
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Figure 2. Sketch of Typical Nodal Layout 


convenience, the network is imagined to be "complete" , even though the shape 

of the material may dictate that some of the mesh boxes are empty. Both 

nodes (boxes or centers) and mesh corners are numbered in a row and column 

system which will be labeled as an m-n system, where m . denotes a row and 

n denotes a column. The m-n mesh scheme may be oriented arbitrarily with 

* 

respect to the physical r-Z coordinate scheme, so that in general the row- 
column description might seem to be merely a descriptive artifice; however, 
in exploiting the simplicities associated with low curvature bodies and heat 
shields it has been assumed in constructing the program that the heated surface 
is at the top of the columns, as indicated in Figure 2. Thus as surface 
recession occurs, the boxes at the top of each column shrink; the other boxes 
remain fixed. This limitation that the heated surface be located at the top 
of each column is merely a convenience procedure exploited for the special 
case of heat shields typical of fairly blunt cones. For high curvature bodies 
such as sharp cones this restriction becomes inconvenient, although Reference 2 
gives an example of the successful use of the code on a very sharp body. 

The side walls of the nodal network are presumed insulated and the back 
wall communicates through a simple heat transfer coefficient law (plus radia- 
tion) to a "reservoir" at ^ Ies - Thus the boundary condition at this face 
does not involve thermochemistry. 

Three other "convenience limitations" have been applied to the layout 

of the nodal grid. First, it is assumed for the purposes of computing thermal 

conductances between nodal centers that the mesh scheme is nearly orthogonal, 

so that conductance may be taken as conductivity times side-area divided by 
* * 

length. Secondly, it is presumed that only one material ablates, and that 
only one material is exposed to the heated surface. Thirdly, it is presumed 
that the principal directions of thermal anisotropy are aligned with the nodal 
mesh. This simplifies computations and reduces input requirements. For those 
applications in which the heated surface intersects the principal directions 
of anisotropy at "difficult" angles this restriction could become a major in- 
convenience and more general schemes would have to be devised for those prob- 
lems . 


That is, the nodal mesh scheme may be above or below the Z-axis line, may be 
oriented in any general direction, and may be "bent" or shaped in any manner 
convenient to the user (subject to the limitations described in the para- 
graphs following below) . 

** 

A suitable more general conductance calculating scheme could remove this 
restriction without any essential change to the program. 
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3.2.2 Geometry 


3. 2. 2.1 Location of Nodal Centers 

The study summarized in Reference 3 showed that to obtain a node 
dropping scheme which produce smooth surface temperature histories, it was 
necessary to avoid associating thermal capacity with the surface temperature. 
In terms of the conventional thermal RC network, this means that a one- 
dimensional scheme would look like Figure 3. All the capacity 



Figure 3. Sketch of One-Dimensional Thermal RC Network 

of the surface node is located in the center of the first nodal zone, and 

half of the thermal resistance of this zone is interposed between the surface 

temperature and the temperature of the first capacity lump. Note that for m 

nodes or boxes the system has m+1 temperature points, in contrast to most 
* 

schemes. Thus for an m x n nodal network scheme, there will be m x n nodal 
temperatures (associated with thermal capacity) plus n surface temperatures 
(not associated with thermal capacity) . 

Denoting the nodal center coordinates as r^ and , the coordinates 
for node m, n with corner coordinates as shown in Figure 4 may be repre- 
sented in terms of the corner coordinates by Equations (1) and (2) : 



Corner 
m+1, n+1 


Corner 
m, n+1 


Figure 4. Sketch of Nodal Center Location For Ablating Material 


* 

The capacity location scheme described here is an improvement over that 
recommended in Reference 3 . 
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c,m+l, n+1 


c,m, n+1 


( 1 ) 


N,m # n 


_ r c,m,n + r c,m+l,n + r 


+ r 




2 ~h Z _l 2 j 2 

_ c, m f n c,m+l,n c, m+1, n+1 0 , 111 , n+1 


(2) 


* 

(Nodes in the ablating material will automatically be centered; the MACABRE 
user has the option to put back-up material nodal centers anywhere in the 
nodal box.) 

3. 2. 2. 2 Path Lengths for Conduction 

In this and the following sections, the nodal center will have a general 

location r XT , Z^ T for illustrative purposes. For computing thermal con- 

N,m,n N,m,n c 

ductances, it will be necessary to have thermal path lengths between nodes. 

Since material properties are associated with each nodal box, it is convenient 

first to consider path segments inside each box. In general there are four 

path segments of interest as shown in Figure 5. The program computes the 

lengths L _ 

^ m,n,B 



Figure 5. Illustration of Path Lengths 

and L ^ in the m direction as the distances between the nodal center 
m. , n , D 

and the centers of the m+1 and m faces of the box. (This is because the 
nodal center is usually on the line joining the centers of these two faces.) 
Thus the paths B and D have the lengths 


* 

Except for the first column, in which nodal points are centered on the first 
side of the box (between corners m,l and m+1,1) under the assumption that 
this line of sides will end at the body stagnation point. 
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r* + IT 

c,m+l,n c, m+1, n+1 _ 


i 2 

+ ^c,m+l,n c, m+1, n+1 _ z 


N,m,n 


N,m, n 


-h 


(3) 


m,n,D 


r . - + r V 

c,m,n+l c,m,n 

2 N,m,nj 


Z _j_ 2 

c,m,n+l c,m,n 


Z N , m, n 


“ft 


(4) 


Similarly, for the n direction path lengths 


m,n, A 


| r c,m,n + r c,m+l,n _ | 

l 2 " N,m,nj 


i z + z 

\ 3 T 

c,m,n c,m+l,n 

Z N,irt,nj 

2 


(5) 


m,n,C 


r c,m+l#n+l + r c,m,n+l _ 


N, m, n 


3. 2. 2. 3 Side areas 


. [ Z c, m+1, n+1 + Z c,m,n+1 ^ \ j 

1 Z N,ra, n) J 


( 6 ) 


Areas of the sides of each box are also required for thermal conductance 
calculations- For the sides lettered as shown in Figure 6, the areas are 


m+1 



m+1, n+1 


Figure 6. Nomenclature of Nodal Sides 
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given by elementary geometry (First Theorem of Pappus) 

A m,n # A " n * r c,m,n + r c,m+l # n^ r c,m+l, n “ r c,m,n^ 


^ Z c,m+l,n Z c,m, J ^ 


A m,n,B ^^Cjirwn + r c, m,n+l^ ^ r c,m. 


•m r 1 3 

n+1 c,m, n' 


+ (Z _ - z 

x c,m,n+l c 


n J 


( 8 ) 


Only these two areas need to be computed for each box, since the areas on the 
other sides of the nodal box are identical with the areas A and B of adjacent 
nodes . 


3.2. 2. 4 Volume 

The volume of the elemental box is by the Second Theorem of Pappus, 

V = [ tt/ 3 ( Z j r (r - r ^ 

m / n L ] c ' m ' n L c,m,n+l c,m+l,n ; 

+ ^ "j i 7 r*~» m 

c,m,n+l c,m+l,nj c,m+l,n c,m+l,n 


( r - r + r 2 -r 2 

c,m,n c,m+l,n+l c,m,n ^ c , m-i-1, n+lj 

+ Z c,m+l,n+l \_ r c , m+1 , n+1 ^ r c, m+1 , n ~ r c,m,n+l^ 

+ Z i I r n 

c,m,n+l [_ c,m,n+l 

C r « ^ \ 

c,m+l,n+l c,m,n' c,m+l / n+l c,m,n 

3. 2. 2. 5 Geometric Effects of Surface Recession 

For nodes adjacent to the heated surface, side B may move due to 
surface recession (ablation) . This recession reduces the thermal resistance 
between the surface point and the adjacent nodal point, increases transverse 
thermal resistances (as discussed below) , and reduces the thermal capacity 
associated with the nodal center of the nodal box adjacent to the surface. 



r 3 — it 3 

c,m+l,n c / m / n+l 
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The program assumes that the surface recession occurs so as to maintain the 
ratios 



Figure 7. Sketch of Surface Nodal Box Undergoing Recession 


a 

a 


2 

1 




( 10 ) 


as shown in Figure 7. The b line joins the centers of the m+1 and m planes 

* 

and serves to define the location of the surface point S . The moving corners 

m+l,n and m+1 , n+1 are located accordingly, and the areas and volumes computed 

as before, except that for these nodes A w _ / n _. 

c m,n,A m,n-l,C 

The nodal point or center is also presumed to move so that is is 

* 

always in the arithmetic center of the nodal box . Thus the nodal center is 
“squeezed" toward the back wall of the nodal box as recession proceeds. The 
path lengths for conduction are adjusted accordingly. 

3. 2. 2. 6 Surface Shape 

As noted above and illustrated by Figure 7, it is most convenient to 

consider the surface points on the neated surface as being always located on 

* 

the nodal column center line , that is , on the line joining the center points 
of the m+1 and m planes ("parallel" to the heated surface) , where m is 
the row index of a nodal box at the surface and m+1 is the local corner index 
of the heated surface (see Fig. 7) . The location of the surface points is 
all the information needed about the surface for those ablation problems with 
a specified, input surface recession rate as a boundary condition (the various 
ablation problem boundary condition options are described in Section 4 below) , 
since for those problems it is most convenient to specify, as input, reces- 
sion history along the nodal center line. The nodal grid as input thus serves 


* 

Except, as noted above, for the first column of nodal boxes, in which the 
nodal points are centered on the first side (Side A of Figure 6) of the nodal 
box and the surface point S is located at the corner m+1,1 (Figure 7) . 
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to define the various nodal column center lines, and the recession history 
for each column then defines the history of the surface points in a perfectly 
straightforward manner. However, another important heated surface boundary 
condition option involves not input recessions but various energy and chem- 
istry information sufficient to calculate surface mass loss from energy bal- 
ance considerations (Option 1, discussed in Section 4.3 below) . This ablation 
calculation does not produce recession rates directly, of course; instead 
it produces rates of mass loss from the surface. It will prove convenient, 
nevertheless , to adhere to the concept of the surface point which moves along 
the column center line as recession progresses. To define the surface point 
motion with computed mass loss rates determined from this general energy- 
balance- determined option requires knowledge of the angle between the local 
normal to the surface and the column center line, since computed mass loss 
may be translated directly into recession along a normal to the surface. Re- 
cession along the normal may be projected into the nodal column center line 
once the angle between these two lines is known. The direction of the sur- 
face normal may conveniently be determined from the slope of the surface, 
that is the surface shape, and of course the direction of the column center 
line is known. 

It is obviously not safe to use the slope of the heated (top) surface 
of the surface nodal box to obtain the surface slope, since in general it is 
neither possible nor always desirable to lay out a nodal grid which will 
"conform" to the "real" surface for the entire problem history. Therefore 
the MACABRE program includes various special curve-fit subprograms which 
compute a surface shape at each time step during the solution after examining 
the layout of the surface points. Referring to Figure 8, one curve fit sub- 
program program computes the 



Figure 8. Sketch of Surface Points 

slope of the surface, dr/dZ , at surface point n as the average of the slopes 
between surface points n-1 and n, and n and n+1. Thus 


dr 

dZ 


n 


l/ r n - r n-l 
2 l Z n “ Z n-1 


+ 


r n+l 2 r n 
Z n+1 ~ Z n 


( 11 ) 
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Note that n-1, n, and n+1 are surface points, not nodal corner points. 
(Problems with only one nodal column have only one surface point. This 
requires the program to abandon the average slope scheme and to use the nodal 
box heated surface slope as the surface slope.) 

Other curve fit routines are available which involve fitted quadratics 
in successive intervals. The User's Manual gives recommendations for the 
choice of fit routines for different body shapes. 

Special formulas are used for the slopes at the first and second 
surface points, however, in order to harmonize the array of surface slope 
values with some basic assumptions made in the convective transfer coefficient 
routines. Since the first column surface point is presumed to be the stag- 
nation point, the slope here is set to a very large positive value. It is 
assumed that the surface is spherical between the first and second points , 
so that the second point surface slope is given by 


where 


dr = l ~ V 2 

dZ 2 2y 



( 12 ) 


(13) 


With surface slope determined, the surface movement AS computed during 
the time step may be projected onto the nodal box center line, and then the 



Figure 9. Sketch of Surface Geometrical 
Relationships (exaggerated) 
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new nodal volume computed, as indicated in the exaggerated sketch of Figure 9. 
Actual surface movement during a time step is limited to a small fraction of 
the nodal thickness to ensure a good approximation to the conservation of mass. 
Note again that the first column is an exception, since the surface point is 
located on the left corner, not in the center of the top surface. 


3.3 INTERNAL CONDUCTION PARAMETERS 

The thermal resistances between nodes and the thermal capacity of the 
nodes are calculated each time interval from the material properties of the 
node corresponding to its temperature at that time. The material properties 
(density (p) , specific heat (c) , conductivity (k) , and emissivity (e) are in- 
put as table look up functions of temperature. Linear interpolation is em- 
ployed for material property determination at temperatures intermediate to 
those tabulated. Constant thermal contact resistances may be specified be- 
tween any or all nodes. The nodal capacities and resistances are calculated 
as follows : 


“m,n, 0 


D C V 

H m,n,0 m,n,0 m,n 


(14) 


m ' n ' A ' e Vn+l,J 


m ' n t c + „m, n+1 , A + R * 

k Q k _ a m, n, B 
m,n,0 m,n+l,0 


(15) 


L m,n,B + ^m+l,n,D + R * 


m, n, B, 0 A , . k « 0 T k ■ , 0 T m,n,A 

m+l,n,B \ m,n,0 m+l,n,0 J 


(16) 


Figure 10 shows the locations of resistances R _ _ and R „ . 

m,n ,A, 8 m,n ,B , 0 


The 



• m-l,n 


Figure 10. Sketch of Thermal Resistance Nomenclature 
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resistances on the other two sides of the nodal box are calculated when the 

quantities for the adjacent nodes are calculated. For nodes adjacent to 

the surface, however, A _ ^ A ... _ generally (refer to Fig. 6 for area 

m,n,C m,n+l,A 

nomenclature) ; for these nodes 


m, n, C 


in , n+1 , A 


m,n,A A k _ a . k , „ 
m,n,C m,n,0 ^n,n+l,A m,n+l,0 


R 

+ m ’ n » B 
m,n, C 


(17) 


It should be noted that anisotropic thermal conductivity values k and k* 
are used in Equations (15) and (16) , respectively. 


3.4 IN-DEPTH CONDUCTION SOLUTION 

The in-depth conduction solution is the explicit finite difference type 

often employed for transient heat conduction analysis. The temperature of 

node m,n at time 0* (T ,) is obtained by application of the finite differ- 

m,n , u 

ence energy balance and rate equations to the nodal volume. 

Solving for T fll , one obtains 


L va,n, 6 ' 


T 

m+1, n, 9 + 
_ R m, n, B, 0 


T 

m-, n.4-1, 9 + 
R m, n. A, 0 


T m-l,n,0 
R m-l,n,B, 0 


T 

m, n-1, 6 
R m, n-1 , A, 0 


T 

m, n. 


0 


R 


m, n-1. A, 0 


R m,n,B,0 R m,n,A, 0 


m-1, n, B, 0 


A9 


+ T 


m, n , 6 m ' n ' 9 ( 18 ) 


In the program this equation is used to obtain "new" temperatures for 
all nodes except for two nodes in each column adjacent to the heated surface 
and for back wall nodes. The nodes near to the heated surface are linked to 
the surface temperature implicitly and hence a special procedure is used for 
the temperature of these nodes, as described in the next section. 

Back wall nodes include a quantity 


hA (T 

m,n,D res 


T ) 

m,n, 0' 


+ ae, A (T 

bw m, n, D v m, n, 0 


T ) 

res' 


inside the brackets of Equation (18) , and have T , Q formally equal to 

m— x f n f y 

zero. 
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The explicit relation (18) imposes the familiar stability restriction on 
the time step size A0 = 0'- 6. The MACABRE program automatically employs a 

conservative stability equation for interior nodes: 


A0 


m,n 


a + R + r + r 

m,n,A m,n,B m,n-l,A m-l,n,Bj 


(19) 


Normally for stability, the input parameter n is less than unity. Surface 
nodes are not considered for time interval calculation, as will be explained 


below; back wall nodes include the terms A 
denominator. 


m,n,D 


(h/2 + 4 a £, T 


bw m , n , 0 


) in the 


The automatic stability criterion calculation may be suppressed for any 
node if the user is sure that the allowed time step for that node will never 
be the minimum one for the system. This saves some computer time. 

Alternatively, the stability criterion calculation can be suppressed 
entirely. If it is not used, then AS must be specified. 


3.5 TEMPERATURES OF SURFACE NODES AND SURFACE POINTS 

Temperatures of surface points are determined either by assignment 
(Option 2) or by surface energy balances described in the next section (Options 
1 and 3) . The surface energy balance determines the new surface temperature of 
the nth column T* with an implicit iteration technique. Stability consider- 
ations dictate that the first and second node temperatures also be treated 
implicitly, and that any transverse heat conduction link (across columns) 
for surface temperatures must be implicit. This latter requirement is a 
complex one to meet as columns recede; hence the surface temperature points are 
not linked transversely. Figure 11 shows the implicit and explicit heat 
conduction paths for two typical situations. 
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Boundary Layer Edge 



Surface Temperatures 
First Nodes 
Second Nodes 


n-1 


n+1 


Implicit 

Links 


Boundary Layer Edge 


V7 JV rrf t /j/jy ' Surface 


n-1 


n+1 


Explicit 

Links 


Figure 11. Sketch of Implicit and Explicit Temperature Links 

in Finite Difference Solution, for Two Typical Situations 

In general terms, the equation for the new temperature of a surface 
node (not a surface point ) is 


T 


m / n / 0 1 


■t 


T m,n+1.9 - T in,n, 9 T m,n-l.e ~ T m,n.8 

R **■ R 

m,n,A,0 m # n-l,A,0 


+ _m-l, n, 9 1 1 m / n, 0 1 + T s,n, 0 


‘ m- 1 , n , B , 0 




m, n, 0 1 


A0 


C m,n,e^ Tm ' n ' 9 


( 20 ) 

This equation links the new first node temperature T m n q « to tlle 

new temperature of the next node down T m _^ n q i • Another equation also links 
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T QI and T - _ Q1 ; this is the energy balance on the second node down, 
m t n , 0 in— l,n,u 

This has the same form as Equation (18) with T Q and T , ^ 0 replaced 

m,n,0 m-l ,n , 0 

by. T and T n Q rin the conductance terms: 

u m,n , 0 m-l ,n r 0 


m- 1 , n , 6 ' 


^m-l^n+lfQ T m-l,n,0 


1 f n f A f 0 


T m+l,n-l,0 T m-l,n,0 


/Xl“l / A / 0 


T m-2,n,6 T m-l,n,9 

■R 

m- 2 ,n ,B , 0 


m,n , 0 1 m-l f n , 9 1 


R m- 1 , n , B , 0 


A 0 


m-l ,n, 0 


+ T 


m- 1 , n , 0 


( 21 ) 


Equations (20) and (21) between them have three new unknown tempera- 
tures at each time step: T Q , , T QI , and T n Ql (where m here 

s,n,U‘ m , n , 0 m- 1 , n , 0 

denotes the row index of the surface nodal box in Column n) . The two equations 

may be combined to eliminate T , A , , leaving a simple linear link between 

m x f n f 0 

the surface point temperature T s n 0 • and the surface node temperature T m n Q t • 
Formally we have a relation between the two unknowns as 


T m / n,e' f s,n , 0 1 ^ 


a l 


T 

s 


,n,0 


.+ 


( 22 ) 


The surface energy balance has the general form (as described in the next 
section) 


convection and chemical energy terms (T ftl ) 

s f n / u 


+ radiation to wall - radiation out (T Ql ) 

s ,n , 0 


m _ rn 

__ s,n , 6 1 m ,n , 9 1 

A m,n,B ,0 


i,n ,40 g cond 


(23) 


where. the parentheses denote f Actional relationship. The equation is written 
on a unit area basis. Relation (22) may be substituted into Equation (23) 
to give the non-linear surface energy balance Equation for T 


s,n, 


convection and chemical energy terms (T Ql ) 

S / XI f t) 

+ radiation to wall - radiation out (T QI ) 

s ,n , 0 


, + b _ 


a 2 T s ,n , 6 


^cond 


(24) 
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When T o, has been found from this equation, the new surface node temperature 
s t n , v 

T m n 6 ' be foun(3 from Equation (22) . For surface energy balance options, 

the method for finding T fiI is described in Section 4 below- (In Option 2, 

S / n / ° 

T is known immediately by user assignment and Equations (20) and (21) 

then determine T Q , and T , 0 . at once.) 

m , n , o m— l , n , 0 
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SECTION 4 


IV-V COUPLING TO SURFACE THERMOCHEMISTRY SOLUTION 
4.1 GENERAL ASPECTS 

The direct computational link between the subsurface solution and the 

surface state solution is effected through the surface energy balance calcula- 
* 

tion. This linkage was discussed from the point of view of the in-depth 
solution in Section 3 above and the energy equation written down in schematic 
terms as Equations (23) and (24) . We now must discuss this equation in more 
detail. The energy balance is illustrated in Figure 12, 



Figure 12. Representation of Surface Energy Terms 
During Ablation 

where the indicated control volume is fixed to the receding surface. Energy 
fluxes leaving the control volume include conduction into the material, radia- 
tion away from the surface, energy in any flow of condensed phase material such 
as mechanical removal, and gross blowing at the surface. Energy inputs to the 
control volume include radiation in from the boundary layer and enthalpy flux 

due to the char flow rate. The final input in the sketch is denoted -q . 

a 

It includes all diffusive energy fluxes from the gas-phase boundary layer. 

If the in-depth response computation were coupled to an exact boundary- layer 
solution, the term -q^ would be obtained from that solution procedure and is, 
of course, a complex function of the boundary- layer structure. If, on the 
other hand, the in-depth response is being coupled to a simplified boundary- 
layer scheme such as a convective film coefficient model, as will be done here, 
then the term -q assumes the form of a correlation equation. 

The computation of the surface energy balance requires from the in- 
depth solution a relation between the surface temperature and the rate of 
energy conduction into the material, <3 con< 3* As noted in Section 3 above, this 
relation derives naturally from the finite difference energy balance for the 

*For Option 1 calculations. The simpler surface options will be 
discussed in Sections 4.5 and 4.6 below. 
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the node just under the surface and is a simple linear equation q cond = a 2 T w + b 
With this information, the surface energy balance considerations allow deter-* 
mination of the thermochemical erosion rate m c and surface temperature T^. 

Itwill be useful to keep in mind that, from this point of view, the purpose 
of the in-depth solution at any instant is to provide information about 
^cond^ T w^ * T ^ e sur ^ ace energy balance equation may be written as 

~ q a + q rad + *c h c " q rad " ^ v) w h w " VX/* ~ q C ond = 0 

in out i (25) 

where 

w - “c - E \ 

l (26) 


The relation g cond = f(T ) delivered by the in-depth solution. 
Other dependencies of interest are 


h c VV (27) 

g rad = q rad (T w^ 

out out (28) 


For the other terms , 


we may write 


in general 


T , — q , q h , ) m h, = functions of boundary- 

w q a y rad w L> r ; l layer-edge enthalpy, 

in pressure, local bound- 

ary-layer solution, 
laws for conservation 
of chemical elements, 
chemical equilibria 
and/or kinetic relations, 
upstream events, m c 


(29) 


If the boundary layer transport aspects of the problem are modeled 
by a film coefficient scheme, then Equation (25) can be normalized on the 
mass transfer coefficient in the customary manner: 




p e U e C M 


q rad " g rad 


in 


out 


p e u e C M 


+ B'h 
c c 


B'h 


/ p u C„ 
L , K .e e M 


\ " 


^ond 

p e u e C M 


0 (30) 
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and relationship (29) becomes 


-«» 


Pe u e C M 


q rad K 

in . \ tA 

:> e u e C M w ' P e u e C M 


functions of bound- 
ary layer edge en- 
thalpy, pressure, 
laws for conservation 
of chemical elements, 
chemical equilibrium 
and/or kinetic rela- 
tions. 


( 31 ) 


The generation of these relationships is the goal of Phase IV, the 
surface thermochemistry solution, to be discussed below. For the present, it may 
be observed that the energy balance coupling proceeds as follows: an initial 

guess of the dimensionless mass removal rate is obtained in some manner. 

With this , the quantities -q a /P e u e c M ' h w ' "'r^ h Ji //p e u e C M' and T w are 
obtained from the surface thermochemistry solution. The quantities h c and 
q rad out are then formulated using the T^ so obtained. The surface energy 
balance is then computed, the Q cond as a function of T w having been provided by 
the in-depth solution. In general, however, the sum of the terms will not equal 
zero, but some error. An iteration procedure is then used to select succes- 
sively better estimates of which drive the error to zero. Experience shows 

that Newton's procedure, in which the derivative of the error with respect to 
B ^ is used to compute the next guess for B^ gives good results. 


4.2 


COMPUTATIONAL APPROACHES 


4.2.1 Pre-calculated Table 

It is evident that each iteration in the search for a surface energy 
balance, if performed as described above, would require a new surface chemistry 
solution/ generally in the near neighborhood of many such .previous solutions. 

If this state solution is to be obtained from a large subprogram (rather from, 
say, simple correlations), this constant repetition of thermochemistry calcu- 
lations would consume an ' ^acceptable amount of computer time. This suggests 
that a tabular approach in which the surface state solutions are done before- 
hand for preassigned B^ values would offer significant computational 
economies. As an example of this approach, consider one body point. Assume, 
for simplicity, that chemical kinetics are not important; then the independent 
variables in Equation (31) reduce to pressure P, boundary layer recovery 

enthalpy H , and B*. Figure 13 represents the space of these three independent 
r c 

variables. The course of a typical transient solution might be as represented 
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in the sketch. As time proceeds, solutions progress 



Figure 13. Representation of Course of Independent 

Variables in Surface History, One Body Point 

through the space of the table independent variables B^, H r , and P. The 
dots in the picture represent solutions satisfying the surface energy balance 
as time progresses. Surrounding these points are a number of other points 
examined during the iterations to satisfy the surface energy balance. For 
any iteration, the solution procedure finds itself, so to speak, within a 
cube formed by the bracketing tabular values of B^, H r , and P. The 
dependent quantities T^ t -q & /P e u e C M' h w and S m r h ^ p e U e C M have been pre- 
calculated for these tabular points? relevant values of these quantities 
for the current iteration values of B^, , and P can then be formed by 

interpolation inside the cube and the surface energy balance equation cal- 
culated. If the energy balance is not satisfied to some preselected degree 
of accuracy, a new value of can be selected and the process repeated. 

The interpolation feature of the tabular approach drastically reduces 

the number of surface state calculations while at the same time allowing 
* 

sufficient accuracy. The precalculated table approach has the following 
advantages : 

1. In parametric studies, tables once generated are useable 

for many different problems, yielding even greater economy. 


For example, a typical 2000 time step problem with the usual 5 iterations 
per time step would require 10,000 surface state calculations, which 
require 1 to 3 seconds each for machines in the 7094 speed class. A single 
enthalpy table, on the other hand, might involve only about 300 state 
solutions (30 B£ values times 10 pressures) , a factor of 30 improvement. 
Multiple enthalpy tables reduce this advantage, but usually only a few 
enthalpies are required. 
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2. For most problems, it is difficult to specify a priori an 

u adequate" array of independent tabular values B^, f H^, and P. 

An examination of the precalculated surface tables before 
execution of the actual ablation-problem-plus -in- depth -solution 
can reveal if there are any "holes" in the tables in areas 
where energy terms are varying rapidly. Desirable table points 
can be added before the in-depth response run. 

3. The surface tables are frequently of independent interest 
in themselves for judging the ablation effects under various 
conditions . 

4. Finally, without the tabular approach, the occasional noncon- 
vergent surface chemistry solution would stop the entire in- 
depth solution process. With the precalculated table approach, 
such solutions are automatically weeded out of the tables 
without damage to the subsequent in-depth solutions. 

On the other side of the ledger, some big disadvantages in the 
precalculated table approach are evident however: 

1. Figure 13, which is a realistic schematic view of a typical 
calculation history, indicates that most of the laboriously 
calculated and assembled surface state points in the table 
are never used in the course of a given solution. 

2. The "mechanical" linkage between the surface state solution 
and the in-depth solution, i.e., the transfer of punched card 
surface state output to input of the in-depth program, leads 
to computing delays and occasionally to gross input blunders 
(wrong decks, missing parts of decks, etc.). 

3. Too much storage is required by the large precalculated table. 

4. The system cannot readily be generalized to accomodate more 
parameters. Both storage requirements and computation time 
grow out of reasonable bounds. 

This last difficulty would prove to be a great stumbling block in 
the ablation code described here, which considers chemical kinetics in 
addition to the other independent variables , H r , and P described above. 
"Kinetics" contribute a fourth parameter in the following manner: The 
chemical state routine considers a number of simultaneous equations including 
element balances, equilibrium relations, and kinetically controlled reactions 
(References 4 and 5) . Each kinetically controlled reaction equation has one 
forward rate coefficient, which in turn comprises one pre-exponential factor 
B' m (for the m^such reaction) and a temperature dependent quantity. All 
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equations are normalized on the transfer coefficient this yields the 

desired independent variable B' in the element balances and also produces 
normalized kinetic rate coefficients B m /P e u e c M - Thus any forward rate coef- 
ficients used in the state routine are input in the form B ra /P e u e c M ' and the 

whole set of B /p u values characterizes the role of kinetics. in a given 
m e e M 

problem the set of B 1 s is fixed; thus really there is only one parameter scaling 
the kinetic effects. Exactly what this parameter is called is not important? 
let it be termed B j_/P e u e c ^ where represents the first of the set of kin- 
etically controlled reactions considered. 

Thus the four problem variables, which would be the independent 
variables in the prepared table, are: 

1. Ablation rate B* 

c 

2. Kinetics parameter B ]/P e u e c M 

3. Recovery enthalpy 

4. Pressure 

Four independent variables are plainly very troublesome. First, 
the savings in computation time are seriously compromised, since the number 
of preprepared solutions has gone up by a large number. Second, machine 
storage requirements are now inconvenient; several dependent variables must 
be stored at each point in an array which might include 30 x 10 x 10 x 10 
= 30000 points. These difficulties in the pre-calculated table approach 

lead to the invention of a revised tabular approach described in the next 
sub-section. 

4.2.2 Direct-Coupled Table 
4. 2. 2.1 General Description 

It is obvious that the difficulties involved in the precalculated 

table approach to a "four-dimensional table problem" could be removed with 

a "close-coupling" procedure which only called for calculation of tabular 

* 

values as they are needed during the actual in-depth solution . At any 
instant, the surface-energy-balance-plus-in-depth response solution procedure 
first examines the corners of the "square" (or "cube", depending on the 
number of independent variables) it finds itself in, to see if the necessary 
dependent quantities have been computed and stored. If not, the in-depth 
routine returns to a master control routine and asks that any missing 
values be supplied. The master control routine determines which corners 
_ 

Another way out, not considered here, is to reduce the number of parameters 
from four to two by ignoring kinetics and by assuming equal diffusion co- 
efficients and C M = C„. For this "classical" simple case, as is well known, 
the edge state does not enter the determination of the surface state, and 
one can use the same two parameter surface table in which and P are the 
independent quantities, for all surface locations and all times. 
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need filling , calls the surface state program to fill them, and then returns 
to the in-depth routine, which continues with the problem solution. 

Such an approach minimizes the number of state points calculated, 
since the four parameters in the list above shrink to only three because 
the second, third, and fourth variables ^ B i/P e u e c j 4 an ^ state) can now 

be directly associated with time as the solution progresses , although when 
this is done axial stations must be distinguished from each other, since 
each station will have different time histories .There is a net decrease of 
one in the number of independent variables. If the tables were to be pre- 
calculated, the three variables would, of course, require separate arrays 
specified in advance, since the true histories of P e u e c M an< ^ edge state are 
consequences of the solution due to body shape effects . With direct 
coupling, B^/p u C M tiie e ^9 e state are effectively known functions of 

time for each station since the body shape is to be calculated at each 
instant. Thus the independent variable list becomes: 

1. Ablation rate B* 

c 

2. Axial station 

3. Time 

Thus the direct coupled approach brings us back to a three-parameter 
system and the same general economics as before. Another computational 
economy is introduced since in general this array will not be completely 
filled in the course of a transient problem. 

Additional economies of computer storage are realized since only 
two table values of time are involved in the calculation at any one instant. 
Instead of having to store a whole box of solutions of the sort shown in 
Figure 13, the program needs to store at any one time only adjacent "sheets" 
of one box as indicated in Figure 14. As soon as the time passes out of the 
interval between the two sheets, the earlier (0 n ) of these two sheets is 
discarded and a new sheet at ® n+ 2 ( tempor ari ly empty of any dependent 
solution values) added. 

Thus the direct coupled approach puts the "more- than- three-parameter- 
table" problem within reach economically speaking, and renders it more 



0 dimension includes 
effects of transient, 
shape dependent edge 
state and P e u e c M 



Figure 14. Sketch Indicating Solution Space and 

Solution "Sheets" Stored at Any Instant 
During Solution 


economical than the precalculated three-parameter table approach. Further- 
more, it reduces computer storage considerably and reduces mechanical handling 
problems of punched tables. 

However, in this direct-coupled scheme two important procedural 
matters have to be considered. These will be described in the following 
sub-section. 


4. 2. 2. 2 Some Dangers in the Direct-Coupled Table Approach 

The general attractiveness of the direct-coupled table approach is 
somewhat compromised by two potential problems: 

1. The user will probably want to specify in advance the array of 

tabular independent values. But what if he makes a poor 

choice, so that the in-depth solution has trouble finding a 
surface-energy balance with the tabular values provided? Can 
a scheme be devised to rescue the problem solution by auto- 
matically inserting new tabular B^ values in unforeseen critical 
or difficult areas of the table? 

2. What can be done about the occasional nonconvergent surface 

state solution obtained from the surface state program? Experience 
shows that it is not possible to guarantee a converged solution 
from complex surface chemistry routines; in any surface table 
array usually some of the points prove to be unobtainable without 
some human intervention in the form of better first guesses for 
the surface state program. In the coupled program, the in-depth 
scheme must be able to adjust to this fact and continue calculating 
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even in the occasional absence of a converged entry in the surface 
state table. 

The code described here does not contain an automatic answer to the 
first problem. A somewhat similar coupled code (Reference 6) used for 
charring materials of very complex composition has shown that such a feature 
is less necessary than might be supposed a priori . In any event, a generally 
worthwhile scheme would be very difficult to devise. The code does contain, 
however, a very valuable feature which successfully meets most of the dif- 
ficulties which arise under point number 1. above. Carbon ablation in air 
provides a worthwhile example. Figure 15 shows the dependence on temperature 
of carbon ablation rates in air. Over a wide range of temperature is for 
all practical purposes independent of temperature; in this range the only 
useful independent variable in the surface thermochemical solution process 
is temperature, since it would be nearly impossible to select a sequence of 
values distributed along the plateau. 

Therefore, the description presented above, in which B^ was used as 
one of the independent variables, should be extended to include temperature 
as an independent variable in place of B^ , at least in some cases. Figure 15 
also shows, however, that temperature is not a good independent variable in 
some instances . 



Figure 15. Ablation Rate B£ Versus Temperature 
for Carbon in Air 

At high temperatures B^ becomes strongly dependent on T , so that is the 
preferred independent variable. 

This general problem is rather typical, it turns out. Consequently 
the single independent variable represented by in Figure 14 has, in the 
code, a more complex dual nature. The variable begins with a string of 
ascending surface temperature (not B^) values (chosen by the user) and 
concludes with a string of ascending B^ values (also chosen by the user) . 


29 


During the solution process , the MACABRE code calls for surface state solutions 
at the minimum (first) value in each string and notes the temperature 
discovered for this point. Temperature points in the temperature string above 
this temperature are ignored. In solving for the surface energy balance at 
a given point at a given time, the code compares the previously formed sur- 
face temperature to the two temperatures at the first B^ values for each 
of the two time sheets currently stored. If is less than either of these 
two values, the surface energy balance iteration process will refer only to 
the user- assigned temperature part of the table (with B^ * s and other dependent 
information at each tabular point discovered as necessary by calls to the 
chemistry routines) . Otherwise, the iteration process will refer to the 
user- assigned B* part of the table (with T w 1 s and other dependent information 
discovered in a similar manner as needed) . 

This dual nature mechanism solves many of the difficulties which night 
otherwise arise under point 1. above. The user must still exercise some 
caution, however. For example, in the case illustrated by Figure 15, the 
code user might begin his sequence of independent variables with temperature 
entries starting at the lowest temperature of interest (say, 530 °R) and 
extending up to (as a minimum objective) the point where B^ begins to rise 
above the plateau value for the lowest pressure to be encountered. It would 
be conservative and preferable to add temperatures even beyond the "rise 
point" at the highest pressure in case these points should be needed during 
the solution. 

After these temperature values, B^’s may be added to span the expected 
range of B^ » s . The user must be careful to pick a minimum B^ slightly above 
the plateau value. A B^ value near the start of the plateau or off the 
plateau in the low temperature region is extremely undesirable for two reasons: 

1. The lowest assigned B^ determines the break between assigned 
temperature solutions and assigned B^ solutions; a low B^ will 
in effect "disqualify" all of the assigned temper ".ture points 
across the plateau and leave a large "hole" in the array of 
surface thermochemical solutions. 

2. For equilibrium calculations, assigned B^ solutions at B^ values 
below the plateau have a high probability of non- convergence 
difficulties in the chemistry routines. 

The second practical difficulty cited above for the direct coupled 
table approach concerns possible non-convergent solutions in the chemistry 
routines. Strictly speaking, even one such solution should serve to halt 
the entire solution process, but fortunately the situation usually does not 
demand such drastic action. The chemistry routines allow a certain number 
of iterations to find the state solution; failure to converge within the pre- 
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set. allowed number of iterations is raire and usually means that the solution 
either is oscillating in the very near neighborhood of the true solution 
or has arrived near to the solution and is proceeding very slowly toward it. 
Therefore , usually the MACABRE code may safely accept a certain number of 
non- converged answers.- Beyond this limit , the computation is stopped as a 
safety measure, since repeated failures to converge are usually due to large 
input blunders. 

4.2. 2.3 The Other Independent Variables in the Direct Coupled Table 

Sub-section 4.2. 2.2 above has described in some detail the user input 
of the T string of independent variable values. There are, of course, two 
additional independent variables, but fortunately these are quite simple. 

The first of these is simply body station or surface point. These 
are kept separate, each point having its own string of T^/B^ values. Thus 
the histories of all the surface points would appear in Figure 14 as a 
number of strings of solution points, each wandering in its own horizontal 
plane ( the time -B^ plane) . 

The last dimension to be considered, the time dimension, is represented 
by a user selected string of time values. The same array will be used for all 
stations. It is not difficult for the user to select a set of values appro- 
priate to the transient nature of a given problem. 

Thus , in evaluating the surface energy balance at any given time , the 
surface energy balance routine refers to stored tabular values of T w and 
enthalpy quantities at the two tabular values of time surrounding the current 
time. If the requisite information has not yet been calculated and stored, 
the FILL routine calls the thermo chemical state program to supply this infor- 
mation. If a "previous time" point (at 8 n in the nomenclature of Section 
4. 2. 2.1 and Figure 14) is being filled in, the thermo chemical state is 

naturally evaluated at stored values of the local independent variables in the 

* 

chemistry solution, pressure and recovery enthalpy, at that previous time. 

If a "future time" table point is being filled in, the FILL routine is obliged 
to obtain the state solution at future local values of pressure and recovery 
enthalpy. If the user is not supplying these values as functions of time 
(and thereby suppressing some of the automatic coupling features of the code; 
see Section 7 below) , then these future values must be estimated. Future 
local pressure is estimated as the future stagnation pressure times the past 
ratio of local pressure to stagnation pressure (i.e., based on the correct 
future stagnation pressure, but the past body shape) . Similarly, future 
local recovery enthalpy is estimated as the future stagnation enthalpy times 
the ratio of the past local recovery enthalpy to the past stagnation enthalpy. 

*The transfer coefficient p u C is also an independent variable, as discussed 
above; if kinetics are being e considered it is handled differently, as will 
be explained below. 



(This approximation refers only to "future time" communication with the state 
routine, for which an approximate future value of H r is quite adequate since 
H r has only a slight effect on the surface thermo chemical state evaluation; 
in the surface energy balance equation evaluation itself, the "properly" 
computed current recovery enthalpy is used.) 


4.2.3 Summary 

Three procedures have been considered for coupling a surface state 
solution to the two-dimensional in-depth solution: direct calculation at 

every iteration, precalculated tables, and direct coupled tables filled in 
as the solution progresses. The third was shown to be the only reasonable 
procedure for a two-dimensional in-depth solution. A suitable coupling 
procedure requires a logic routine to determine which table points need 
filling during the surface energy solution process , plus some means of 
specifying the values of the table independent variables. In the code, the 
user specifies a priori the sequence of B* values be considered appropriate 
for each body station, plus a series of time points at which solutions are 
to be performed. 


4.3 THE CONVECTIVE SURFACE ENERGY BALANCE EQUATION 

Up to this point, we have been able to discuss the IV-V link between 
the in-depth solution and the surface thermochemical state solution, which 
is effected through the convective surface energy balance operation, without 
discussing the surface energy equation itself except in its general nature, 
given by Equations (23) , (24) , (25) , (30) , and (31) . It is now appropriate 
to consider this equation in more detail. As noted previously, the surface 
linkages are built upon a transfer coefficient approach: transfer coefficient 

equations are used for mass diffusion in the element balance equations and 
a transfer coefficient type of surface energy balance equation is employed 
in the IV-V linkage. 

The proper choice of energy equation constitutes a subject far too 
large and difficult to discuss completely in this report. The MACABRE code 
uses one of three different forms of the convective energy equation, 
depending on user choice. In ascending order of complexity, and descending 
order of "well foundedness" , these equations are 


(1) For equal mass diffusion coefficients for boundary layer 
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°w g rad in 


molecules and for = C H (Le = 1) : 


H - ( 1+B 1 ) h + B'h + 
r w C C p el i e C H 


Foe 4 

p e U e C H w 


^cond 

p e u e C H 


(32) 


32 



This equation is a "standard" (see, for example. Reference 7) one 
and requires no discussion. 


(2) For equal diffusion but C' ^ C„ (Le ^ 1) : 
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(33) 


This equation supposedly accounts for ^ C H effects , and doubtless 

does so if the boundary layer is frozen. The accuracy of Equation (33) 
for chemically active boundary layers has been much discussed and the question 
remains open (References 7, 8, 9) . 

(3) For unequal diffusion: 
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(34) 


This equation, introduced in Reference 10, appears to be useful for 
modeling unequal <3. f fusion effects, at least in frozen boundary layers, and 
currently is being studied in detail to establish its validity (Reference 9) . 

All three of these equations (32) , (33) , and (34) serve to provide 

expressions for the diffusion energy flux -q Q in Equation (30) . The con- 
stituents of -q , which include h. , 2 Z* hT w , and 2 Z*i h^ w 

a edge x e 1 ^ 1 

gas 

are generated by the chemistry routines along with h w and T . The subroutine 
FILL, which calls the chemistry routines, then combines the quantities in 
the appropriate fashion as dictated by Equation (34) or one of the simpler 
forms (33) and (2). For example, suppose the FILL routine discovers a 
requirement for a state solution at a tabular and for a given pressure 
and recovery enthalpy. FILL calls the state routines, controlled by EQUIL , 
which perform, the state solution, evaluate key property values, and return 
to FILL values T , h , £Z* h . Tw , h , 2 Z* hT w . FILL then, referring 

W W 1 1 ” 1 1 
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to the temperature T , looks up the ablating material enthalpy h c and con- 
structs the term 


Y 


edge 

gas 


M 


Z i 

(Z* - Z* )h. w + B'h - B'l 

1 I X c c 


and stores it as one of the dependent quantities of interest. The temperature 
is stored as the second dependent quantity of interest. (On the assigned 
lower part of the table, EQUIL returns and this is stored.) 


4.4 CONVECTIVE SURFACE ENERGY BALANCE ITERATION PROCEDURE 

At every surface point at which a convective energy balance is to be 
obtained, the energy balance subroutine must 

(a) decide whether T^ or B' is the proper independent variable 
to iterate 

(b) select a trial value of this independent variable 

(c) refer to stored dependent values of energy quantity Y 

and temperature (or B^ if in first part of table), inter- 
polating as necessary on B^ and 6, and filling in needed 
table points as necessary by calling FILL 

(d) look up « w (T w ) and (T w ) , q rad (6), local F, plus C M , 

in 

C^, and ( 0) (which in many cases are generated by special 

subroutines described in Sections 8 and 9 below) 

(e) formulate the left-hand side of Equation (34) , note departure 
from zero (error e) 

(f) if error is too large, correct trial value of independent 
variable, go to (c) - 

Generally this process converges very fast, usually within three 
of four iterations. Step (f) , the correction step, is based upon the 
Newton- Raphson procedure: at the same time the error e is being computed, 

the derivative of the error with respect to the independent variable is 
also formed, using derivatives of the tabular quantities Y and T^ (or B^) . 
The correction in the independent variable is then computed directly: 


AB^ = -e/ (de/dB^) 

or 

AT w = -e/(de/dT w ) 


(35) 

(36) 
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4.5 


AVOIDING THE SURFACE STATE SOLUTION 


The convective boundary condition described so far in this section 

does not include all problems of interest. Consequently the surface boundary 

condition treatment has been formulated as a triple option. The first of 

these (called Option 1) constitutes the convective boundary condition described 

in Section 4.1-4. 4 above. A second option (Option 2) merely refers to input 

time-dependent quantities for user assigned values of surface temperature 

and recession rate. A third option (Option 3) retains the energy balance 

aspects of Option 1, but eliminates the convective heating terms and does 

not allow surface recession. This option is useful for "cool-down” or 

"soak -out" portions of problems, and can also be used for assigned heat 

flux problems in general (through the q , . variable) if surface recession 

2 T sq in 

is negligible. The heating option may be changed at the user's choice at 
any time for a given column or surface station, and the heating option assign- 
ments may vary in any manner around the body at a given time. 



SECTION 5 


SURFACE STATE CALCULATION - PHASE IV 


5 * 1 INTRODUCTION 

Section 4 above has described the convective surface energy balance 
boundary condition and described the use made in this boundary condition 
calculation of computer subroutines for evaluating various heated surface 
thermochemical state variables, chiefly enthalpy and molecular composition. 
For a restricted range of problems it would be possible to avoid the use 
of chemistry routines for this calculation and instead rely on enthalpy 
correlations and the like, but for the present coupled code it was strongly 
desired to have a general ablation program which could treat an Y non-charring 
ablating material exposed to any arbitrary environment. Fortunately , 
suitable thermochemistry routines already existed which could be coupled to 
the in-depth response code; the program reported here did not develop these 
chemistry codes. Since these codes are fully described elsewhere (Refer- 
ences 4 and 5) this section will seek only to give a general overview of the 
physical problems treated by the chemistry routines. 

The thermochemical routines as a package are referred to as ACE, 
after Aerotherm Chemical Equilibrium code. This is an inappropriate 
acronym since (1) the routines consider non-equilibrium chemistry as well 
as equilibrium, and (2) in the coupled code no routine has the name ACE; 
the controlling routine is called EQUIL. For historical reasons, however, 
we shall continue to call the thermochemistry package ACE. 

5.2 USES OF ACE IN MACABRE 

The ACE package is called upon for a number of computational tasks 
in the coupled code. Its main use, of course, is in generating surface 
thermochemical state solutions. As will be discussed below, these are open 
system calculations involving interactions between edge gas and the main 
ablating material. The surface energy equations (33) and (34) include 

terms like h and EK. h^ w or ZZ* hT w , however , which indicates another 

w edge gas X e 1 ^e 1 

needed thermo chemical calculation: the molecular make-up of the boundary 

layer edge gas must be identified and certain properties evaluated (such 
as Z| ). Each call to ' ACE for a surface state solution in fact involves the 
following comput ations : 
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(1) 


a closed system equilibrium calculation on edge gas 

alone to determine K. and Z* ; this calculation is 
1 _ 1 _ £ 

done at assigned pressure and recovery enthalpy . 

(2) an open system ablation calculation at assigned pressure, 

p u (for kinetic effects) , and either B* or T 
e e M c w 

(depending on location in independent variable array) ? 
this calculation may include kinetic effects 

(3) a frozen evaluation of h and ZZ* h. w at 

edge x e 1 

gas 

temperature T^. 

The quantities h , Zz* hT w , h , and ZZ* hT w are obtained 

w edge 1 e 1 w V 1 

gas 

from ACE and employed as described in Sections 4,3 and 4,4 above, 

(The ACE routine may also be used by the pressure distribution 
and transfer coefficient routines to generate certain thermodynamic proper- 
ties and transport properties of the boundary layer gas. These properties 
are stored for future use in a Mollier Chart (as will be described in more 
detail in Sections 8 and 9 below) in which P and h g are the independent 
variables; the required ACE calculations are therefore closed system equili- 
brium calculations of the edge gas at assigned pressure and enthalpy.) 

The following sections describe the manner in which calculations 
are done in the ACE package. 

5.3 PHYSICAL PROBLEMS TREATED BY THE ACE PACKAGE 

The ACE package of subroutines treats a wide range of thermo- 
chemical problems. Since these problems are complex, the capabilities of 
the ACE package can perhaps best be illustrated by a set of descriptions 
of increasing generality. The following sections describe such a graded 
series of problems. 

5.3.1 Closed System in Equilibrium 

Closed systems are defined as those for which the relative 
amounts of chemical elements are prespecified. The edge gas state cal- 
culations are of this type. 

Consider K chemical elements, , introduced into a previously 


Good arguments can be made for using the static enthalpy instead, but this 
proves inconvenient in the code. The difference in results is small in 
any case. 



evacuated container. In general, these elements will interact to form a 

* 

number of chemical species , (gas phase) and (condensed phases) . 

If enough time has elapsed so that thermodynamic and chemical equilibrium 

is established, the thermodynamic state of the system, including the relative 

amounts of chemical species present, is completely determined if two inde- 

* * 

pendent thermodynamic variables are known . This condition may be stated 
mathematically by examining the governing equations for such a system, 
and showing that the number of independent equations is equal to the number 
of unknown quantities . 

Relations expressing the formation of the gaseous chemical species 
from the gaseous chemical elements may be written as follows : 


2 C ki N k ' N i (37 

K=1 

Similarly, formation of condensed phase species from the gaseous elements 
is written: 

K 

2 c ja N k + N * (38 

K=1 

In the above, represents the number of atoms of element k in a molecule 

of species i (gas) or species £ (condensed) . 

If the gas phase species are assumed to individually behave as 
thermally perfect gases, then the equilibrium relation corresponding to 
reaction (37) is 


P. 


l 



k=l 


K . (T) 
Pi 


In P. 

l 


E 

V=1 


c ki ln p k 


In K . (T) 
Pi 


(39) 


where P k denotes partial pressure and K pi (T) is the equilibrium constant 
for the formation reaction (37) of species 1SK . For each candidate condensed 
— * ~ 

"Chemical species" as used here includes molecular, atomic, ionic, and 
electron species. 

■ * 

Duhem 1 s Theorem . 
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phase species 


K 

-E 


k=l 


C]U lnP k = m V (T) 


( 40 ) 


where 

= indicates the existence of the condensed phase species 
in equilibrium with gas phase species, and 

< indicates that the condensed phase species will not 
be present in equilibrium. 


For each chemical element introduced into the system, the conservation 
of atoms dictates that the amount of any element k in the gas and con- 
densed phases (regardless of molecular configuration) must sum to the total 
amount of element k in the system. Mathematically, this may be written, 
for each element k, as 


Mass fraction 
of element k 
input to the 
system 


w E + w E v* 

i=i *=i 


( 41 ) 


where M is a composite system molecular weight* defined by 
I ^ L 




■ E t * E «' nsg.lfTa"' ) 

i = l £=1 


( 42 ) 


and where X p is a mole fraction of condensed phase species l defined as 


_ molecules of condensed species & 
x Jl “ total gas phase molecules x 


( 43 ) 


* 

This is the molecular weight appropriate to the ideal gas equation of state 
if condensed phases are present. 
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In addition, for the gas phase species, there exists the requirement that the 
partial pressures must sum to the total system pressure 


I 


E 


p i 


= P 


(44) 


Mixture 
related to the 


thermodynamic properties, 
species concentrations by 


such as specific enthalpy, 
equations of the form 


are 


h = 



p i h i 


+ 


1 


y x h 

1=1 


(45) 


Consider now the number of independent equations for the system. The 
number of gas phase equilibrium relations (39) is equal to the number of 
gas phase species I minus the number of elements K (because Equations (39) 
are trivial when i=k) . In addition, there exists a relation such as (40) 
for each of the L candidate condensed phase species in the system. Note 
that the system temperature is contained implicitly in Equations (39) and (40) 
through the temperature dependence of the equilibrium constants. There are 
K conservation of elements Equations (41) , one for each atomic element 
introduced into the system. The requirement that the partial pressures sum 
to the system pressure (44) contributes one additional equation. For any 
additional thermodynamic properties of the mixture (enthalpy, entropy, etc.), 
there exists equations such as (45). 

Consider next the variables appropriate to this formulation of the 
problem. The relative concentrations of the I species in the gas phase are 
given by the P^'s and the amounts of the L candidate condensed phase 
species are given by (most or all of which may be zero) . In this formu- 
lation, the composite system molecular weight, 5 h is also a variable. There 
are one each of the mixture thermodynamic variables T, P, h, s, etc. The 
number of variables and available independent equations may be summarized as 
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NO. OF SUCH 

EQUATION 

NO. OF SUCH 

VARIABLES 

VARIABLES 

NUMBER 

EQUATIONS 

p i 

I 

(39) 

I - K 

x * 

L 

(40) 

L 


1 

(41) 

K 

p 

1 

(44) 

1 

T 

1 

of the type 


h,s, p, • • - 

n 

(45) 

n 

total 


total 


variables 

I+L+n+3 

equations 

I+L+n+1 


Thus, there are two less equations than there are variables, and so if two 
independent variables are specified (e.g., P and T) in addition to the 
elemental composition, then closure is obtained and the chemical and thermo- 
dynamic state of the system may, in principal, be determined. 

The ACE program performs this determination. That is, to determine 

the equilibrium thermodynamic and chemical state of any closed system, one 

needs only to furnish the ACE program with the elemental composition, the 

* 

candidate gaseous and condensed phase species, and two thermodynamic pro- 
perties. One of these properties must be pressure, and the other may be 
either temperature, enthalpy, or entropy. Given this information, the ACE 
program will calculate and output the mole fractions of each candidate 
species, the temperature, enthalpy, entropy, density, effective molecular 
weight, equilibrium and frozen specific heats, the isentropic exponent, and 
a few other quantities of potential interest. 

Ac already noted, MACABRE only calls upon ACE for assigned pressure 
and enthalpy solutions. The edge state is saved momentarily while the 
open system calculation is performed. 


The ACE program must be supplied certain basic thermodynamic data (described 
in References 2 and 5) for each candidate species to be considered in a given 
system. This data is contained in a three card set, one set for each species. 
A certain amount of judgement is required on the part of the user relative 
to which candidate species should be included in a given system. Frequently 
this judgement is avoided by simply inputting data for all species containing 
combinations of the input elements. 



5.3.2 Open Systems in Equilibrium - Simplified Case 

The basic theory underlying the treatment of open systems may best be 
illustrated by examining the equations expressing the conservation of chemical 
elements and energy at the ablating surface. If the boundary layer is charac- 
terized by equal diffusion coefficients, and if no material is removed from 
the surface in a condensed phase (i.e., no mechanical erosion or liquid layer 
removal) , then these equations take on a particularly simple form for equi- 
librium systems. This simplified situation will be considered in this section 
and these considerations will be extended to more generalized cases in 
Sections 5.3.3 to 5.3.5. 

Consider first the fluxes of chemical elements (k) entering and 
leaving a control surface affixed to the ablating surface. The solid material 
may be visualized as moving into this surface at a rate s. If it is assumed 
that no material is being removed in a condensed phase, then the surface 
and the fluxes of the k th chemical element may be illustrated as 
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surface 
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Figure 16. Ablating Surface Mass Balance 

Terms superscripted by a tilde U) represent the total mass fraction or flux 
of element k, independent of molecular configuration. Thus 



i=l 


a k K i 
1 w 


(46) 
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( 47 ) 



where k pertains to element k, i pertains to species i, and < is the 
mass fraction of element k in species i. Fluxes of element k away from 
the surface consist of boundary layer diffusion and gross motion of the fluid 
adjacent to the surface due to the injection flux m c# 

From the above sketch, it is seen that conservation of chemical 
elements requires that 


jfc w + (pv> A, = 


m c K k 


(48) 


Summing Equation (48) over all elements k yields the total mass continuity 
equation (for the case when there is no condensed phase material removal) 


( Pv ) 


W 


(49) 


An important fundamental of the present mathematical modeling of the ablation 
process is the expression of the diffusive heat and mass fluxes in terms of 
a transfer coefficient formulation. 

In this formulation, the diffusional term -j, is written 

K w 


-3 


, = p u C (K 1 

k e e m k^ 

w e 




(50) 


Utilizing this transfer coefficient formulation (50) for the diffusion 
flux in the elemental balance Equation (48) yields 


P u C m (K. 
*e e M k 


■ v 


<» v >» K k 


m K, 
c k. 


(51) 
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In terms of dimensionless ablation rate and solving (51) for the total 
mass fraction of element k at the wall yields (for equal diffusion coefficients 
and no condensed phase material removal) 



B i K k + \ 

c e 


1 + B* 


(52) 


Given the relative amounts of chemical elements specified by (52) , the 
chemical and thermodynamic state of the gases adjacent to the ablating surface 
may be calculated from equilibrium relations. 

Since the gases are in equilibrium with the ablating surface 


- E 

k=l 




in P, = in 


V (T) 


(53) 


if i represents the surface species, and 

K 

" Z c k* in p k < ln V (T) 


k=l 


(54) 


for all other candidate condensed phase species. In the present formulation, 
(54) may be thought of as being used to identify surface species i for 
which (53) applies. The equilibrium relations for gas phase species is of 
the form 


in - 


K 


Z 

V=1 


in P k = in K ± (T) 


(55) 


and the partial pressures must obey the relation 


I 


z 

i =i 


P. 

i 


= P 


(56) 


The elemental mass fractions adjacent to the surface, , 
related to the species partial pressures, P., by relations 


of (52) are 
such as 



- \y cl . p . 

i=l 


( 57 ) 
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Thus, if P is known and is specified (this may be varied parametrically 
as will be discussed subsequently) , and if T w and P i are unknowns, then the 
number of unknowns and equations available may be summarized as 


UNKNOWNS 

NO. OF SUCH . 
UNKNOWNS 

EQUATION 

NO. 

NO. OF SUCH 
EQUATIONS AVAILABLE 

p i 

I 

(55) 

I-K 

Nr 

K 

(53) 

1 

ft 

1 

(52) 

K 

T 

w 

1 

(56) 

1 



(57) 

K 

Total 

Unknowns 

I+K+2 

Total 

Equations 

I+K+2 


Thus, closure is obtained and, in principal, the temperature and 
chemical composition of the gases adjacent to the surface may be determined 
if P and are specified. From the pressure, temperature, and chemical 
composition, the calculation of other thermodyanmic properties (enthalpy, 
etc.) is straightforward. 

The ACE package of programs has the capability to determine the 
chemical and thermodynamic state of the gases adjacent to an ablating sur- 
face in a fashion similar to that discussed here. In the coupled MACABRE 
code, the ACE group is typically called with pressure P known and 
assigned as an independent variable by the user. The temperature T^ and 
other state data are then determined. Alternatively, T^ may be assigned 
(in the lower temperature parts of the tables, as discussed in Section 

4. 2.2. 2 above) and B' is determined. 

c 

5.3.3 Open Systems in Equilibrium - Unequal Diffusion Coefficients 

The discussion of the previous section was limited to open systems 
with equal species diffusion coefficients and no removal of surface material 
in the condensed phase. These simplifications were made in order to render 
the basic theory easier to explain. While this simple model is reasonably 
accurate for many ablation situations, these assumptions are inappropriate 
for others. For this reason, calculations performed by the ACE package 
(and hence the MACABRE program) are not restricted to any of these simplifi- 
cations. The extension to unequal diffusion coefficients will be discussed 
first, in this section. 
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The unequal diffusion coefficient model has already been briefly intro- 
duced in Section 4.3 above, in which Equations (33) and (34) may be compared. 

In simplistic terms, the model replaces the mass fraction during forces K. 

* . . e 

with modified driving forces . A similar modification must be made m 

the ACE treatment. The diffusive flux of element k at the wall is now 

modeled by 



p e u e C M 




(58) 


In (58) , Z* is, in effect, a weighted average of the mole and mass fractions 
of element k. The Z* are defined by 


K 



i=l 


(59) 


Z* 
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z I K i" Y 


zYk?“ y 

/ i x 


i=l 


(where y « 2/3, see 
Ref. 10 ) 


(60) 


V F i 
E V F i 

i=l 


( 61 ) 


where the F^ are diffusion factors defined by the following relation for 
the binary diffusion coefficients 





(62) 


where D is a constant for a given temperature and pressure and the depend 
weakly on temperature. The JK j must obey (6 2) in order for the boundary 
layer species diffusion equations to reduce to a form from which (58) can 
be inferred by similarity arguments. Reference 11 demonstrates that the 
binary diffusion coefficients for a variety of chemical systems are accurately 
correlated by (62) . This reference also shows that a reasonably good cor- 
relation equation for the F^ is 
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where 


ffiref “ 23-4 


and e « 0.431 


( 63 ) 


when D is taken as the self-diffusion coefficient of 0 2 . Additional dis- 
cussion relative to this unequal diffusion coefficient formulation is con- 
tained in Appendix A, 

Consideration of unequal diffusion coefficients affects the surface 
elemental balance relationships which are needed to. determine the equilibrium 
thermochemical state at the surface. Substituting (58) into (48) yields 


».w z ; w - kj * - vv * v, 


(64) 


and the "unknowns" here are K^. and Z* , each of which may be expressed 
in terms of the species partial pressures 

I 

2^ i i 

i=l \ 

s»r L, Vi 

E Vi 1-1 


K k = 


i=l 


(65) 


and 


E SciVl 


Z,* 


i=l 


« I 


* - ^ 


Z \ p i /p i 


~ 3c Z c ki p i /p i 


9 7= 


i=l 


i=l 


(66) 


where ^ g is the mean molecular weight of the gas phase and F is a mean fT 
defined as 


Z¥i 

i=l 

I 

Z^v F i 
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Substituting these expressions into (64) and utilizing the definition of 
yields an expression for the species partial pressures at the surface in 
terms of quantities at the boundary layer edge and in the material 


i=l 


P iw + P 


i=l 


C ki P ^/ pY = 


1W 


P»L - 

' V lz » 




( 68 ) 


Note that (68) reduces to (52) when the diffusion coefficients are equal. 

When performing unequal diffusion coefficient open system calculations, 
the ACE program utilizes (68) rather than (52) as the elemental mass balance 
equations. Other than this, the solution philosophy is essentially as 
discussed in Section 5.3.2. The diffusion factors utilized in the ACE pro- 
gram may be calculated in three ways, at the user's option 

a. diffusion factors F^ may be input individually for 
each species i 

b. diffusion factors may be calculated according to (63) 
with the user specifying the reference molecular 
weight , ?/( ref # and the exponent e 

c. if the user does nothing special, the program will 
calculate F^ according to (63) with ^ ref - 23.4 and 
e = 0.431. 


The actual program input for these alternatives is discussed in Reference 2. 
It should also be pointed out that the diffusion factors have an effect on 
the other transport properties calculated by the ACE program, and these are 
briefly discussed in Appendix A. 

For unequal diffusion coefficients, the transfer coefficient formula- 
tion for the surface energy equation is Equation (34). Consistent with (34), 
for unequal diffusion coefficient problems, the surface thermochemistry 
quantities prepared by the ACE program include the quantity 


r 




w 


in addition to the quantities previously discussed. Note again, that for 
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equal diffusion coefficients 


t - 1 V'" - “v 

i=l W i=l W 


Similarly, the quantities 


and 


i=l e 


w frozen 
edge gas 


are computed after reference to the previous automatically computed and 
stored edge solution. 


5.3.4 Ope nSy s terns In Equilibrium - With Condensed Phase Surface 

Material Removal 

Considerations to this point have been based on the assumption that 
all the mass leaving an ablating surface is in the gas phase and this flux 
has been denoted by (pv)^. However, for many ablation situations of interest, 
some of the material leaving the surface is in a condensed phase - e.g. , 
liquid layer removal or "mechanical ablation". The condensed phase removal 
affects both the surface mass and energy balances (Reference 10) . For the 
surface elemental balance, there is an additional flux term 


Z "Si, 

1=1 

leaving the surface. This results in the surface elemental balance (in 
terms of the species partial pressures - analogous to (68) having the form 


'LSci 

i=l 
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p e U e C M 


(69) 
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for each atomic element k. Summing over all elements yields the obvious net 
mass continuity relation 

L 

* c = (pV) w + XA 

A=1 (70) 

Similarly, there is an additional energy flux term 

L 

S.= l 

leaving the surface. The surface energy balance equation becomes (analogous 
to (34) 


p e u e C H (H r h w* frozen + C M 
edge gas 


E T 

(Z* - Z* )h. w + B’ h 

1 11 c c 

e w 


E- 


m 0 H 0 - (B 1 ) H 
l i w w 


+ q , . - a eT 4 = q , 

^rad in w ^cond 


(71) 


The condensed phase surface material removal model currently incor- 
porated in the ACE program is termed a "fail temperature model". According 
to this model, a "fail temperature", T f a j_i_£/ may be pre-assigned to any 
candidate condensed phase species &. This species will then be removed 
from the surface (in a condensed phase) at a rate m^ if the surface tempera- 
ture is greater than, or equal to, the fail temperature for this species. 
Additional comments relative to the computational aspects will follow a 
brief discussion of the fail temperature concept. 


The fail temperature is the temperature at which the material "fails", 
i.e., the temperature above which the material will not remain affixed to 
the surface. For melting solids for which a liquid layer is unstable, the 
fail temperature is simply the melt temperature. For thermochemically 
eroding solids , the fail temperature might be set at some temperature 
characteristi c of the material yield temperature at the appropriate external 
loading. Also, the ACE program has a provision for inputting a maximum 
fail temperature for all condensed species. Physically, this maximum fail 
temperature might correspond to the fail temperature of the substrate. 

When failing occurs, an additional equilibrium relation is avail- 
able for each condensed phase species l for which m^ > 0 . Also the surface 
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temperature must be consistent with the assignment fail temperatures: 

T - T.- . - for the surface species 
w fail 

for any species for which > 0 

These additional equilibrium relations, surface temperature restrictions, 

and assigned fail temperatures are sufficient to solve for the additional 

unknowns m^ introduced into elemental mass balance relations. The enthalpy, 

, of the condensed phase being removed is determined from the species 

thermochemical data. If T f < T melt , H £ is taken as the enthalpy of the 

solid, and if T, - T 14 _, H 0 is taken as the enthalpy of the liquid. The 
f ai l me it & 

convention for associating fail temperatures with each or all candidate 

* 

condensed phase species will be specified in Reference 2 . 

In communicating back to MACABRE, the ACE package amalgamates the 

terms (pv) H and 

w w L 

^2 V 1 *, 

1=1 


by suitably modifying the returned value of H . Thus FILL constructs the 
quantity Y as described in Section 4.3 without specific reference to failing 
events . 


5.3.5 Non-Equilibrium Effects in Open System Calculation 

To calculate the equilibrium state of a chemical system, information 
relative to all chemical reactions is not needed. This fact permits a 
significant simplification in the problem formulation, and those features 
of the ACE package discussed thus far take advantage of these simplifica- 
tions. For some systems of interest, however, the effects of reaction kin- 
etics may not be neglected. A general solution of complex problems for which 
reaction kinetics effects are important is potentially difficult for at 
least two reasons: a) there are significant computational and bookkeeping 

problems associated with the analytical treatment of mixed equilibrium and 
nonequilibrium systems, and b) for many systems of engineering interest, 
the rate controlling reactions are not well known and/or rate constants 
for these reactions are unavailable. The ACE package of routines effectively 

*For carbon ablation in air, the fail temperature concept is of almost no 
interest. Graphites loaded with metals, however, tend to form thin 
oxidation resistant oxide coatings which generally ablate by melting. 

The ACE package will properly identify any such oxides formed, and the 
fail temperature mechanism allows them to be ablated realistically. Fail 
temperatures are useful for metal heat shields, of course. 
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surmounts the first of these problems as discussed in Reference 4 and as 
will be briefly summarized subsequently . However, difficulties falling into 
the second category above frequently preempt any such "exact" analytical 
treatment. In these cases, various approximate treatments can often yield 
useful information. These approximate treatments essentially consist of 
not allowing certain species, or groups of species, to react with an other- 
wise equilibrated system. A few examples of this sort of approximate 
treatment will be discussed in the following paragraphs. These discussions 
will be followed by a brief summary of the general kinetics option of the 
ACE program. 

5. 3.5.1 Removal of Selected Species from the System 

The simplest, and perhaps the crudest, way to approximate certain 
kinetic effects is to simply remove from the equilibrium system those 
species whose formation is, in reality, suppressed by reaction kinetic 
effects. Computationally, these species are removed merely by removing 
the thermo chemical data cards for that molecule from the species thermo- 
chemical data deck. (Reference 2 describes this deck in detail.) 

5. 3. 5. 2 Isolated Species 

Certain chemical species may be frozen individually at any given 
concentration. This treatment is sometimes useful for species which are 
relatively nonreactive because of chemical kinetic effects. For example, 
consider the flow of a C-O-H gas over a relatively low temperature ablating 
carbon surface. Isolation of at its boundary layer edge concentration 

would provide a more realistic surface thermochemistry table than if full 
equilibrium were assumed (since the water-gas reaction at the surface is 
relatively slow at lower temperatures) . In the ACE program, this isolation 
of species is achieved by treating that particular species as an atomic 
element, an element with a fictitious atomic number so that it cannot 
react chemically with any other members of the chemical systems. Com- 
putationally, this requires specification of the input of the appropriate 
quantity of this "element" (accompanied by appropriate reductions in the 
amounts of other actual elements) , and the alteration of the species 
thermochemistry data card set for that species so that it contains only 
one "atom" of the fictitious "element." 
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5. 3. 5. 3 Isolated Subgroups of Elements 

in some practical applications, there exists a subgroup of elements 
which are in equilibrium with each other, but the subgroup is not in chemical 
equilibrium with the system as a whole. The isolation of a subgroup of 
elements can be achieved by considering these elements to be different "ele- 
ments" with properties the same as, but with atomic numbers different from, 
the corresponding actual elements (e.g. , add 100 to the atomic numbers of 
the elements in the isolated s ub group) . For this type of calculation, the 
subgroup to be isolated is considered to be one component, and the com- 
position of this component is specified as the actual relative elemental 
composition of the subgroup to be isolated, but in terms of the altered 
elements. In order that these altered elements may equilibrate with each 
other, the species thermochemical data card set must contain a set of species 
which may be formed from these altered elements, i.e., species with the same 
atomic number convention. 

5. 3. 5. 4 Low Temperature Surface Equilibrium Suppression 

As previously discussed, equilibrium is an increasingly poor assump- 
tion at low surface temperatures (e.g., below 2000°R) in some ablation 
problems. However, in most surface thermochemistry tables, lower tempera- 
ture ranges are usually only of minor interest (e.g. , during a very short 
initial period of a transient ablation problem). Nontheless, the surface 
tables must frequently extend to these lower temperatures in order to be 
compatible with the heat conduction and ablation solution. For these 
cases, the ACE package contains an option by which surface (heterogeneous) 
equilibrium will be suppressed and surface recession nulled (m c = 0) below 
some input limit temperature for solutions in the assigned temperature 
part of the surface tables. The gases adjacent to the surface are assumed 
to be in equilibrium , but heterogeneous equilibrium is not required for 
this solution. For many cases, surface thermochemistry tables generated 
in this fashion are closer to reality than full equilibrium solutions down 
to low temperatures. The input convention for implementing this option 
is discussed in Reference 2. 

5. 3. 5. 5 Full Treatment of Kinetics 

The previous discussion has described various approximate treatments 
of kinetics effects which are often useful when a general solution, including 
an "exact" treatment of reaction rate effects, is either impossible or 
impractical. The ACE program does, however, possess the capability of 
accurately treating a number of kinetically controlled reactions in otherwise 
equilibrated systems (either closed or open) . The remainder of this section 
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will be devoted to a brief description of this "exact" treatment of kinetics. 

The inclusion of kinetically controlled chemical reactions is 
accomplished by removing equilibrium relations (39) from the set of equa- 
tions for certain species participating in reactions that are to be 
kinetically controlled. These equations are replaced by kinetic rate 
equations (of the Arrhenius form) for each kinetically controlled reaction. 
This is accomplished by first identifying the primary reactive species in 
the reactions which are to be kinetically controlled, and second, by allow- 
ing these species to be created or destroyed only via the kinetic rate 
equations. This approach requires a relabeling of species to be considered 
in the kinetically controlled reactions. These species are called pseudo- 
elements since they behave like elements except that they may be created 
or destroyed at rates specified by the reaction rate equations. 

Computationally, the inclusion of kinetics results in the addition 
of a rate-of-creation or destruction term to the elemental balance equations 
for these pseudo-elements. This adds additional unknowns tc the system 
equal in number to the number of species whose concentrations are kineti- 
cally controlled, i.e., the pseudo-elements. The relative creation and 
destruction rates of all pseudo-elements in a given reaction are related 
by stoichiometry, however, so the number of additional unknowns remaining 
is equal to the number of kinetically controlled reactions. The reaction 
rates, from which the pseudo-element creation or destruction rates derive, 
are given by 




TT P v j* _ _L_ TT P *> 

11 j j K pm j 3 


moles/ft 2 sec 


(72) 


for each kinetically controlled reaction m of the form 


y P R N. + 

Z-T D m D 3 rn j 

j j 


(73) 

where the sums and products are over the species and pseudo-elements , 

U are the stoichiometric coefficients, and superscripts R and P denote 
reactants and products respectively. In (72) , K is the equilibrium 
constant for reaction (73), and k Fm is the forward rate constant. In the 
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present formulation, k_ is represented by an Arrhenius type function 

i? m 


$ 

k_ = BT m exp (-E /KT) 

Fm m * w (74) 

where E m is the activation energy, <J> m is the temperature exponent, and 
is a constant representing a multitude of phenomena. E m , 4> m , and 
must be input to the program for each kinetically controlled reaction. 

These constants should be based on experimental data. The uncertainty in, 
or unavailability of, these data for many chemical systems of interest in 
ablation problems frequently represent a significant constraint on the 
application of the kinetics option of ACE program to these problems. More 
detailed discussion of the computational treatment of kinetics is presented 
in References 4 and 11. Reference 2 describes the input details for the 
kinetically controlled reaction data. 

5.4 SOME CODING DETAILS 

The computational procedures are employed in the ACE program to solve 
the equations set forth in Section 5.3, briefly discussed here. Considerably 
greater detail is presented in Reference 4 and in particular, Table I of 
that reference. 


5.4.1 Basic Solution Technique 

The basic solution technique may be illustrated by considering, for 
example, an open system with unequal diffusion coefficients, no condensed 
phase material removal, and no kinetics (i.e., as discussed in Section 
5.3.3) . For this system, the basic equations defining the problem are the 
elemental conservation equations (67) , the total pressure equation (56) , 
the reaction equilibrium equations (55) , and one heterogeneous vapor pressure 
relation (53). The table of Section 5.3.3 shows that there are as many 
knowns as unknowns in these equations so closure is obtained. 


Summarizing these equations, 

B *E c ki p i + F XXiV F i - ? d 
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(75) 
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P = 0 


(76) 
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( 77 ) 


(78) 


The number of unknowns could immediately be reduced by I-K through 
the direct substitution of (77) , as solved for P^ , into the other relations. 
It proves advantageous, however, to continue to treat the full set of 
equations , and to subsequently utilize this substitution during the itera- 
tive convergence procedure. The solution of these simultaneous nonlinear 
algebraic equations is based on Newton-Raphson iteration. Since this pro- 
cedure is accelerated by casting the equations into a more linear form 
(via transformations, substitutions, etc.) it is well to examine the set 
of equations above. With the boundary layer edge, char and pyrolysis gas 
composition given as well as the B', (75) and (76) are linear relations 

between the and ^ providing that F is reasonably constant. In con- 
trast, (77) and (78) are linear relations between the in P^, £n P^ and 

in K . , the latter variable being approximately linear in 1/T. 

P^* 

The ACE program takes advantage of this situation by treating those 

species which are significant in the mass and pressure balances (75) and 

(76) in terms of P^ and the less significant species in terms of their 

in P. . 
i 

The Newton-Raphson procedure, as applied by the ACE program, can be 
summarized by considering a set of equations of the general form 


f j ^ X 1 ' x 2 7 * • ' ' x i ' * * • ) — ^ 


At any point in the solution procedure there exists a set of estimates, 
x* , for all the variables which will in general not satisfy all of the 
relations and will lead to a non-zero value of the f ^ , namely, e ^ . The 
Newton-Raphson method proceeds to "drive" these errors toward zero by 
evaluating the change in each unknown variable, Ax^ , which would reduce 
all the errors to zero if the functions, f ^ , were linear. The linear 
approximation is based on the current values of the unknown variables 
and the corresponding array of values of the partial derivatives 3fj/3x^, 
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Thus 


df . 


i-L 


9f . 


3^ dx i 


which is locally correct and is integrated to 


(Af.)* (AX i )] 


in the linear approximation. The solution of (80) is 


dx 




(79) 
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( 81 ) 


where the array of partial derivatives appearing in (81) is simply the matrix 

inverse of the array in (80) . In the ACE program the formulation of the 

f 

partial derivatives uses the variables, Jin P^, Jin T and In 7ft and (81) 
yields, for example. 
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which if taken as linear all the way to solution yields 


Jin 



(A Jin P^ * 


^/9 Jin P. 
“ 2 ^ 


(-Ej)* 


(83) 


since the desired change in the functions is simply the negative of the 
error. An equally exact relation obtained from (82) is 


dP. = 

i 



3 in P i 


af. 


(84) 


^The choice of Jin P. permits a matrix reduction by the use of the simple 
algebraic substitution, previously mentioned after (78) , prior to matrix 
inversion. 


57 



which if taken as linear all the way to solution yields 


P** - 


p i 


= AP* 




( 85 ) 


The ACE program uses (85) for all species significant in mass balances and 
(83) for the others. 

5.4.2 Restriction on Corrections 

The set of correction (Ax.)* can be thought of as a vector in the 
space of the independent variables which is added to the current vector 
approximation x£ to yield a new estimate x|* . Experience has shown that 
it is frequently unwise to proceed along this correction vector the full 
amount indicated by (83) or (85) . Rather, it is better to proceed a limited 
way, although preserving the same direction. At other times, it is expedient 
to depart from this vector, and seek another based on freezing the value of 
some variable and eliminating a corresponding equation. 

The scaling of the correction vector is such as to limit changes in 
the partial pressures of major species to increases of one order of magni- 
tude and decreases of three orders of magnitude, and changes of temperature 
to approximately 20 percent. 

Molecular weight, temperature and condensed phase concentration 
corrections are frozen and a new correction vector generated if the initial 
set of corrections indicate excessive temperature or molecular weight 
excursions, a contradictory temperature change, or negative corrections on 
newly introduced condensed species. 

The formulation of these and other scaling and freezing criteria 
is an essential feature of the ACE program. Because of these features , 
convergence is virtually assured for well formulated, physically unique 
problems . 

5.4.3 Base Species 

The discussion of Section 5.3 described the equilibrium reaction 
equations as equations giving the formation of a species from atomic elements . 
Thus the reactants are elements and the products are usually molecules . 

This scheme has the advantage of formal simplicity, since the stoichiometric 
coefficients needed in the equilibrium equations are given directly by the 
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product species chemical formula. This scheme can have computational dis- 
advantages, however, since the atomic elements are frequently not present 
to any great extent in the equilibrium system. This in itself results in 
no disadvantage. If, however, a molecule (e.g. , CO) dominates more than 
one mass balance (e.g. , C and 0) t loss of significant figures can slow or 
defeat convergence. 

It is more desirable to write the equilibrium reactions (37) and 
(38) as well as the mass balances in terms of reactant species which are 
actually present in appreciable amounts. These species are termed "base 
species" (from Reference 12) since all other species are taken to be formed 
from them. 

The ACE program selects the base species from the candidate species 
thermo chemical data which are input as specified in Reference 2. The 
program automatically selects as base species the first set of species 
satisfying the requirement that (1) all other species may be formed from 
this base species set and (2) that no balanced reaction can be written 
involving only base species. One base species may be considered to represent 
each element. Thus, the base species are established by the order in which 
the user inputs the candidate species thermo chemical data. The calculation 
of the stoichiometric coefficients and equilibrium constants appropriate to 
any set of base species is handled automatically by the program. 

While the selection of base species is not usually critical, certain 
computational economies may be realized if the base species are those species 
which occur in relatively high concentrations in the system. Also the treat- 
ment of a chemical kinetics (Section 5.3.5) may have an influence on the 
selection of the base species. Therefore the Users 1 Guide (Reference 2) 
recommends certain input procedures for forcing the selection of desirable 
base species . 
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SECTION 6 


COUPLING THE TRANSFER COEFFICIENT CALCULATION TO THE 
SURFACE THERMOCHEMISTRY SOLUTION (III - IV COUPLING) AND 
TO THE IN-DEPTH TRANSIENT SOLUTION (III - V COUPLING) 

Section 5 above makes it plain that the transfer coefficient approach 
to modeling diffusional fluxes is built into the surface thermochemistry 
solution process, since the chemistry routines use as the measure of mass 
injection and the concept depends naturally on the transfer coefficient 
idea. An additional, more explicit, link to the transfer coefficient model is 
the scaling of kinetic effects on P e u e c M - The pre-exponential factors B^ 
are input directly to the coupled code, and these are normalized on P e ^ e C M 
inside the ACE package, as noted in Section 5 above. This requires that a 
p e U e C M value communicated directly to the ACE package with every call 
for use in normalizing and B^'s. Other coupling links appear in the surface 
temperature dependence and the "blowing reduction" of the transfer coefficient 
these are discussed in Section 9 below. 

The III - V coupling between the transfer coefficient and the in- 
depth solution occurs rof course, in the transfer coefficient equation for 
the surface energy events. The relevant equations have already been pre- 
sented in Section 4.3 and require no further comment. 
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SECTION 7 


COUPLING THE FLOW-FIELD (OR BOUNDARY LAYER EDGE STATE) SOLUTION 
(PHASE II) TO THE SURFACE STATE SOLUTION (II - IV COUPLING) 

AND TO THE TRANSFER COEFFICIENT SOLUTION (II - III COUPLING) 

Before discussing the Phase II and III calculations themselves , it 
will be of interest to examine the coupling links between the Phase II 
solution, which refers to the specification of the boundary layer outer 
edge state around the body (pressure, static enthalpy, recovery enthalpy, 
and transport properties) , and the Phase III and Phase IV solutions. 

The boundary layer edge state influences the surface thermochemistry 
state solution, as has already been discussed in Section 5 above, through 
the direct communication of local pressure and recovery enthalpy to the ACE 
package. The edge state in turn is influenced by the current body shape, 
as will be discussed in Section 8 below. 

The boundary layer edge state influences the transfer coefficient 
calculation through the distributions of pressure and enthalpy around the 
body, as might be expected. This is discussed in Section 9 below. 

All of these coupling links are of interest, but it is also useful 
at times to be able to suppress some or all of them. Such coupling suppres- 
sion is particularly desirable in parametric studies attempting to study 
one effect at a time, rather than coupled effects, and in studies in which 
it is desired to match experimental data in some respect (such as measured 
surface temperatures) . Therefore, the MACABRE code has been written to 
allow varying degrees of decoupling. All of the following options exist 
for each surface point on a body at anytime: 

a. Surface temperature and recession rate discovered by 
energy balance, or 

b. Surface temperature and recession rate assigned as 
functions of time by user (Option 2) 

Under (a) we have 

a. 1 . Full surface energy balance with convection and thermo- 
chemistry (Option 1) 

a. 2. No convection, no recession, assigned heat flux (Option 3) 


Under (a.l.) we have 

a. 1.1 Local pressure and/or local recovery enthalpy assigned by user 
as a function of time 

a. 1.2 Local pressure and/or local recovery enthalpy computed by 
appropriate subroutines automatically 

For both (a. 1.1) and (a. 1.2) we have 

a.l.(i).l Local heat transfer coefficient assigned by user as a 
function of time 

a.l. (“) .2 Local heat transfer coefficient computed automatically 
by appropriate subroutines 

For both (a.l.(j).l) and (a.l.(-~).2) we have any degree of "blowing reduction" 
correction to the transfer coefficient. 
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SECTION 8 


THE BOUNDARY LAYER EDGE STATE SOLUTION 


A complete and accurate analysis of the boundary layer edge state is 
a very complex problem requiring simultaneous consideration of the body shape, 
shock shape, the shock layer flow, and the boundary layer flow. For many 
engineering purposes, however, a detailed analysis of this sort is not pos- 
sible, and certain simplifying assumptions are necessary which bring the 
analysis to a more tractable level of sophistication. Since the application 
for the analytical technique being discussed here is primarily re-entry 
vehicle nosetips and heat shields, the first simplifying assumption is that 
the entire region of interest is immersed in a boundary layer edge flow which 
has passed through the detached bow shock at or near the stagnation point. 

Thus, it is not necessary to calculate a shock shape and shock layer flow 
since the shock layer gas of interest is assumed to pass through a normal 
shock. The next simplifying assumption is that the boundary layer edge gas 
state is determined by an isentropic expansion from the body stagnation 
condition. With the edge gas entropy fixed at the stagnation value, the 
complete edge gas thermodynamic state is determined by the specification of 
one more property, typically the local static pressure. Thus, the deter- 
mination of the edge gas state reduces to a problem of first calculating the 
stagnation entropy, which is relatively straightforward, then determining 
the pressure distribution around the arbitrary body shape, and finally 
evaluating the other edge state properties of interest by carrying out the 
isentropic expansion to the known local pressure. The techniques used to 
carry out these steps are described in the subsections below. 

8.1 THERMOCHEMICAL COMPUTATIONS 

Thermochemical computations are required in a number of instances in 
the coupled code, including boundary layer edge state calculations. As 
mentioned in Section 5, these calculations are carried out by the ACE routine. 
For boundary layer gas properties , it has been found most convenient to store 
ACE solutions for the requisite thermodynamic and transport properties in 
tabular form and supply this table of properties to the program as input. 

The option also exists of merely specifying the desired values of the inde- 
pendent variables in this "Mollier Chart' and calling for ACE calculations to 
determine the required properties. The properties at the gas state of interest 
during the course of a solution are found by interpolation between tabular 
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values. The precalculated Mollier Chart approach has a number of advantages 
similar to the precalculated surface tables: reuseability , external evalua- 

tion of the table accuracy and sufficiency, and elimination of non-converged 
solutions. However, care must be taken to guarantee that interpolation within 
the chart will yield sufficiently accurate state properties. 

The Mollier Chart uses pressure and enthalpy as independent variables. 
Dependent variables found at each (P-h) combination are temperature, density, 
Prandtl number, viscosity, entropy, and isentropic exponent. Linear interpo- 
lation is used everywhere within the chart, however pressures and densities 
are converted to logarithmic form before this linear interpolation is carried 
out. Since the computer code dimensions limit the number of pressure-enthalpy 
combinations which can be input, it is recommended that the user check the 
accuracy of interpolated solutions before going ahead with problem solving. 
This can be accomplished by plotting the properties of interest, determining 
the region of most likely error within the linear interpolation approximation, 
and comparing an interpolated point with an "exact" solution from the ACE 
program. 

The edge state solution procedure is typically started by using the 
specified values of stagnation pressure and stagnation enthalpy to determine 
the entropy level from the Mollier Chart. The static pressure around the 
body must then be found, as discussed in the next section. 

8.2 STATIC PRESSURE DISTRIBUTION 

The technique used to calculate the static pressure distribution 
includes a variation of Newtonian theory in the subsonic flow near the stag- 
nation point, matched to a Prandtl Meyer expansion downstream of the sonic 
point. The methods used in these regions are explained below. 

8.2.1 Subsonic Region 

The pressure distribution in the subsonic region is found by a varia- 
tion of the technique suggested in Reference 13 , in which the Newtonian 
pressure distribution relation is modified to include blunt bodies with sonic 
corners. The basic Newtonian pressure expression is 


P 


+ (1 - PJ cos 2 6 


( 86 ) 


where barred quantities are normalized by the stagnation pressure, and 6 is 
the angle between the body axis of symmetry and a normal to the body surface 
at the point of interest. The sketch below will help to clarify this nomen- 
clature . 
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Figure 17. Diagram of Nosetip Showing Nomenclature 

A class of nose tips of particular interest are smoothly contoured 
blunt bodies having local slopes such that the flow is subsonic everywhere 
ahead of a convex sharp corner and the flow expands to the supersonic con- 
dition at the corner. For these shapes, the Newtonian pressure distribution 
is grossly inaccurate since it depends only on local slope. It is obvious 
that the pressure should vary smoothly from the stagnation value P = 1.0 
at the stagnation point to the sonic point pressure P = P* at the sharp 
comer. However, for a sphere-cone shape as shown in Figure 17 above, the 
Newtonian description of pressure would call for a constant P determined by 
the surface slope along the entire cone surface. It was found in Reference 13 
that a modification of the Newtonian distribution involving an empirical 
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expression for pressure on a flat disc of the same diameter (sonic point 
diameter) would greatly improve pressure predictions. The suggested pressure 
distribution expression, modified for nonzero ambient pressure, is 



(87) 


where P™ is the flat disc pressure distribution expression, s is the surface 
running length, is the nose radius, and R max is the maximum of either the 
nose radius or R* , as defined in Figure 17. It is necessary to specify 
expressions for the sonic point pressure P* and the flat disc pressure 
in order to use .this expression. It is in this last expression that the 
present analysis differs from that of Reference 13 . 

The sonic point pressure has been chosen to be described by the ideal 
gas expression 

( 88 ) 

The isentropic exponent y is evaluated at the stagnation point and is assumed 
to remain constant over the region of interest. The assumption of a constant 
y introduces a negligibly small error into the analysis. Another quantity 
of interest is the location of the sonic point, s*. For the sharp cornered 
bodies discussed above, s* is merely the distance to the corner. However, as 
these sharp cornered bodies ablate or for other smoothly-contoured bodies , 
the sonic point must be located in order to use Equation (87 ) . The Newtonian 
pressure distribution is employed here to give the surface angle at the sonic 
point of 



(3* 



(89) 


Thus the surface running length s* can be found by using the known surface 
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slopes and running lengths along the surface and interpolating between them. 
The flat disc pressure distribution used in Reference 13 is 


P FD = 1 - e- X (l-P*) 
X = 4 Vj>n(^-) 


(90) 


It will be shown that the heat transfer coefficient at the stagnation point 
depends upon the stagnation point velocity gradient (see Equation (120)). 
This gradient in turn depends upon the second derivative of pressure with 
distance according to the following relation 


— 1 = 1 

J - — 

d 2 P 

ds | s+o 1 

V p o 

ds 2 


Differentiating Equation (87) twice yields 



s+o 



(91) 


(92) 


which for appLoaching infinity yields a zero velocity gradient at the stag- 
nation point if Equation (90) is employed as it stands. To correct this un- 
realistic result, the data of Reference 14 were evaluated to obtain a flat 
disc stagnation point velocity gradient expression 


du 

ds FD , s-*o 



(93) 


By combining Equations (91) , (92) , and (93) , the following result is 

obtained 



s+o 



(94) 
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This expression provides a more correct stagnation point pressure distribution 
The flat disc pressure distribution, Equation (90) , has been modified to blend 
smoothly with the new second derivative requirements, 

P FD = 1 - e" X (1-P*) - ± [(f*) 2 - e' X ] (95) 
and the definition of X has been changed slightly: 


X = 5 V«.n(-*) (96) 

s 

Equations (95) and (96) are the final forms of the flat pressure distri- 

bution relations which are used in the present analysis. 

8.2.2 Downstream of the Sonic Point 

Using the expression given in Equation (87) above, the local static 
pressure is exactly equal to the Newtonian pressure for all body shapes at 
the sonic point. The Newtonian description of pressure is then used down- 
stream of the sonic point until the Prandtl-Meyer match point is reached, 
at which point a Prandtl-Meyer expansion expression is used for the remainder 
of the running length. The match point is determined by the requirement of 
continuity of pressure and pressure gradient along the surface: 


dP 

ds 


Newtonian 


dP 

q s 


Prandtl-Meyer 


(97) 


It is more convenient to work in terms of the angle 8 and dimensionless 
pressure : 


dP 

dB 


Newtonian 


dP 

dB 


Prandtl-Meyer 


(98) 


The Newtonian pressure expression can be differentiated to give 


dP 

dB 


Newtonian 


(1-PJ sin2 B 


(99) 
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For a perfect gas, the Prandtl-Meyer pressure gradient is 


dP 

d3 


Prandtl-Meyer 


ypm 2 


while the relation between pressure and Mach number is 


( 100 ) 


P = ^1 + M 2 ^ 1 Y ' (101) 

Equating the two gradient expressions yields the following equation for the 
match point: 


e 


M 


1 

2 


PyM 2 ) 


(1-PJ 



( 102 ) 


This is an implicit equation for 3^, since P and M are functions of 3. The 
solution is found by a Newton-Raphson iteration procedure. 

Pressures on the body downstream of the match point are found by 
carrying out a Prandtl-Meyer expansion from the match point conditions to 
the local surface slope. The governing equations in this region are 
Equation (101) above relating pressure to Mach number, and 


e tan_1 V^T (M2 “ 1) " 

■ i ^^r + 


- 


(103) 


which relates Mach number to local surface angle. Once again, the equations 
are implicit since M is a function of 3 / therefore a Newton-Raphson iteration 
procedure was devised to solve these equations . 
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8.2.3 Examples of the Accuracy of the Method 

The accuracy of the pressure prediction technique described above 
is demonstrated with four sample problems. The first of these is a 60° 
half angle sphere cone with a sharp corner as shown in Figure 18 . With this 
shape, the location of the sonic point is fixed at the sharp corner, there- 
fore 6*, s*, and R max are known from the geometry. The pressure prediction 
for this shape is compared in Figure 18 to the data of Reference 15 , and 
also tc the Newtonian pressure distribution for the same shape. Agreement 
is seen to be excellent over most of the surface, which is to be expected 
since the empirical portions of the theory were invented to solve this very 
problem. The shortcomings of Newtonian theory are also readily apparent in 
this figure. 

Figure 19 demonstrates the utility of this pressure prediction tech- 
nique on a spherical shape. For spheres, the Newtonian description of pres- 
sure in the subsonic region is generally very good, and th_ case is no 
exception. The differences between Newtonian and the currently proposed 
theory on a sphere in the subsonic region are very slight, however, and both 
methods give good agreement with the experimental data. 

Figures 20 and 21 demonstrate the agreement between theory and 
experiment for a blunt and a sharp ellipse, respectively. Again, there is 
very little difference between Newtonian and the current theories, with 
either giving satisfactory agreement with the data. 

In summary, the present theory offers good agreement with Newtonian 
theory for spherical and elliptical shapes, but also provides accurate 
descriptions of pressure over flat-disc- type shapes where Newtonian theory 
fails entirely. This makes the theory particularly suitable for the large 
cone angle blunt shapes being considered for planetary entry. 
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Figure 18. Pressure Distribution on 



60 Sphere Cone with Sonic Corner 
















8.3 


EVALUATION OF OTHER PROPERTIES 


Once the pressure is known at any given body point and the entropy 
is assumed equal to the stagnation value, other gas properties of interest 
may be found by referring to the Mollier Chart tables as described in 
Section 8.1. Some care is necessary in carrying out this operation, however. 
The Mollier Chart lookup subroutines have been written to use pressure and 
enthalpy as the independent variables in the lookup procedure since these 
are the known quantities at the stagnation point. In order to minimize the 
error resulting from linear interpolation, properties at other body points 
where pressure and entropy are known are found iteratively by estimating 
enthalpy, looking in the tables at the known P and estimated h, determining 
entropy s, comparing with the known s, revising the h estimate, and so on. 
This procedure prevents a compounding of errors which could arise if one 
were to interpolate on P and h to find s , then interpolate P and s 
to find h. 
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SECTION 9 


THE TRANSFER COEFFICIENT EVALUATION 


An energy integral technique is used to evaluate the nonablating 
heat transfer coefficient along the surface of an arbitrary body shape which 
starts at a stagnation point. The energy integral technique is completely 
discussed in Reference 17, and will be only briefly summarized here. 

In the energy integral technique, the flat plate relationship be- 
tween energy thickness and heat transfer rate is assumed to be valid for all 
types of flows. This relationship is then utilized to solve the energy 
integral equation for an arbitrary flow. To outline this method, it is use- 
ful to establish a few basic definitions. 


heat transfer coefficient: 


% 


p’WV 


(104) 


energy thickness: 


r 6 

£ w / pu 

r h w J p e u 


H - H 
e 


e e \ H r h w , 


dy 


(105) 


f H - h 

e w 

H - h 

i r w 


(106) 


reference enthalpy: 

h* = 

0.231i 

e 

+ 0.19H r + 

0.58h w 

(laminar) 


h' = 

0.361i 

e 

+ 0.19H + 

r 

0.45h w 

(turbulent) 

(107) 


where 

H r = 

u e 

h e + E -f 


(108) 


r = (Pr’) 1 / 2 (laminar) / 
r = (Pr 1 ) 1 / 3 (turbulent) j (109) 


and all primed quantities are evaluated at 
the reference enthalpy condition. 
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It will be recalled that for flat plate flows with no initial boundary 
layer, the heat transfer coefficient expression is 


/-> 1 — 

a 

C H ” 

(Pr') 2/3 (- 

where 

a = 0.332 ! 


m = 0.5 

or 

a = 0.0296 


m - 0.2 


laminar 


turbulent 


( 110 ) 


Combining this with the definition of C H ' results in 


p e u e < H r"V 


a(p') 1 m (u') m 
p e u™ x m (Pr ' ) 2/3 


( 111 ) 


which relates heat flux to x for a flat plate. For the flat plate with nearly 
isothermal walls, the energy integral equation is 


d<J> q w 

dx p e u e (H r _h w> (112) 


Combining Equations 111 and 112 and assuming small variation in h^ compared 
with H r “h w f the flat plate energy integral equation can be solved to yield a 
relationship between energy thickness and x: 


x 


(Pr 1 ) 


2/3 


O)’ 


-,1/1-m 


(1-m) <f> 


a (p'/p_) 


(113) 


By combining Equations 113 and 111, we have an explicit relationship between 
heat flux and energy thickness: 


P 


eWV 


.(Pr') 


a(p’/p e ) 

2/ 3 


1 

1-m 


(114) 


This expression is now assumed to be valid for planar bodies or bodies of revo- 
lution, with arbitrary pressure gradient. 
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The more general boundry layer energy equation in integral form is 


dx 


rp u dx (rp e u e J + (H -h ,) dx (H r h w } 


p u (H -h ) 
e r w 


(115) 


Combining Equations 115 and 114 , the resulting form ca. e integrated directly 
to give 


1 

1-m 


frp u (H -h )~1 

L e e r w'j. 


rP e u e (H r -h w ) 


1 

1-m 


■ — 1 — I ■ 

[ rp e u e (H r~V] 1 "” x i [ 


arp'u e (H r -h w ) (y ' ) 
Pr 2/3 (p'u ) m (l-m) m 


1-m 

d\ 


(1-m) 


(116) 


where X is a dummy variable representing running length. This relation between 
and x can be combined with Equation 114 above to give the general relation 
between heat flux (or heat transfer coefficient) and x. For laminar flows 
starting from a stagnation point, the final expression is 


p e u e C H 


0.332p'y'ru e (H r -h w ) 


(Pr ') 2/3 j fo r 2 p'y'u e (H r -h w ) 2 dx ! 1/2 


'117 ) 


Similarly, for all turbulent flow, 


P U C TT 

e e H 


0.0296 p’u e (y 'r) 1/4 (H r -h w ) 1/4 
(Pr ' ) 2/3 | y T p , u e (y') 1/4 r 5/4 (H r -h w ) 5/4 dxj 


1/5 


(118) 


Finally, for transition from laminar to turbulent flow at x = , 

x 

/ 

(Pr') 5/6 < 0.9028 


p e u e C H = (0 ' 0129 'u e (y , r) 1/4 (H r - h w ) 1/4 

2 V /8 

p'y'u e (H r - h w ) aX J + 0.0161 


{C - 


4/3 


/ 


r 5/4 , , ,,1/4 

ftf P'u ip ) 

(Pr ') 5/6 e 


(H r " V 5/4dX 


1 1/5 


(119) 
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For a stagnation point, we use the laminar expression and take the limit as 
x + o with 


r = x 



p 1 u 1 = const. 


H-h = const, 
r w 

The heat transfer coefficient reduces to 


p e U e C H = 


(Pr 1 ) 


664 
2?3 


p'y ' 


du 
e 

dx 


1/2 


( 120 ) 


Consistent with the pressure distribution discussed in Section 8.2.1, the 
stagnation point velocity gradient can be expressed as 


which becomes 



( 121 ) 


( 122 ) 


In the present version of the code, only the laminar heat transfer 
expression. Equations (117) and (120), are included. For programming purposes, 
it is most convenient to write the heat transfer expression in the form 


where 


p e u e c t 


0.332 

(Pr') 2 / 3 


( p' P 


.jl/2 



(123) 


u* 

l* 


du , 

e / stagnation 

ds V point 


(124) 
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or 


u* = u„ 
e 

s 

f P ' U ' ( H r " h w ) 2 r 2 u e ds 

£* = - 9 . 

P'y' <H r -h w ) 2 r 2 u e 

The running integral term for l* away from the stagnation point is evaluated 
by treating each linear term in the integrand as if it varied linearly be- 
tween body stations where its value is known. That is, p', y', H r -h w , r, 
and are all approximated as linear functions of s between each pair of 
body stations, thus making direct integration possible. The value of the 
integral at each body station is stored and merely added to fhe increment 
in £* at each new body station. All primed quantities are evaluated by re- 
ferring to the Mollier Chart lookup routine with the known pressure and 
reference enthalpy. 

An example of the accuracy of the energy integral method is provided 
in Figure 22, where results of the present theory applied to a sphere- 
cone configuration are compared with results from Aerotherm's multicomponent, 
nonsimilar, boundary layer analysis program (Reference 18). The two programs 
were supplied with identical pressure and wall temperature distributions over 
the length of interest. It is apparent that the energy integral approach gives 
satisfactory results for the nonablating heat transfer coefficient for at least 
a sphere-cone type of body. 

For use in the MACABRE code, the nonablating heat transfer coeffi- 
cients computed with the analysis outlined above are corrected to account for 
the effects of ablation ("blowing"). If we denote this nonablating coefficient 
(that is, for m = 0) as P e u 0 C H ^, then both data and simple analyses (see, for 
example. Reference 7) indicate that the dependence of P 0 u e c H on m is given 
fairly accurately by 


away from 

stagnation 

point 


(125) 


(126) 



(127) 


where 


cp = 


2 Am 

P e u e C H_ 


( 128 ) 


This correction is built into the surface energy balance operations of the 
MACABRE code. 
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SECTION 10 


OVERVIEW OF COMPUTER PROGRAM EVENTS 


The preceding sections have described the individual computing phases 
of the coupled computer code and the linkages between the phases. This sec- 
tion presents two simple flow charts for the code to clarify further the 
operation of the code. Figure 23 gives an overall chart of code events, and 
Figure 24 shows some details of the surface energy balance sub-package. 
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INCREMENT MAIN ITERATION 



Figure 23. Overall MACABRE Code Flow Chart 
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CALLED FROM MACABRE 

i 

COMPUTE ALLOWABLE TIME INCREMENT AS 
INFLUENCED BY SPACING OF TIME TABLE 
ENTRIES, SPACING OF MANDATORY SUR- 
FACE STATE SOLUTION TIMES, RECESSION 
RATES, OUTPUT TIMES 

♦ 

COMPUT CURRENT SURFACE SHAPE, SLOPES 


MAIN LOOP OVER SURFACE POINTS (COLUMNS) 
SET COLUMN INDEX 1=1 



CONTINUED ON NEXT PAGE 


Figure 24. 


Surface Energy Balance Operations Flow Chart 
(Subroutine SURFB) 




1 

IF HAVE JUST ARRIVED AT NEXT MANDATORY 
SOLUTION TIME, LOOK IN RELEVANT TIME 
TABLE FOR THIS STATION TO FIND FUTURE 
LOCAL PRESSURE, FUTURE LOCAL RECOVERY 
ENTHALPY, AND FUTURE LOCAL TRANSFER 
COEFFICIENT; ALL OF WHICH WILL BE RE- 
QUIRED IN FUTURE-TIME CALLS OF THE 
SURFACE STATE PACKAGE. IF SOME OF 
THESE QUANTITIES ARE NOT ENTERED AT 
THE FUTURE TIME, USE SPECIAL RULES TO 
ESTIMATE THEM FROM CURRENT DATA 

I 

CORRECT TRANSFER COEFFICIENT FOR BLOWING 
REDUCTION 

CHECK LOWEST B' ENTRIES IN SURFACE TABLE 
FOR VOIDS; IF ANY, CALL ACE THERMOCHEM- 
ISTRY PACKAGE TO FILL IN DEPENDENT VARI- 
ABLE DATA INCLUDING T 

w 

i 

DETERMINE IF SURFACE ENERGY BALANCE 
SOLUTION SHOULD BE OBTAINED IN THE 
LOWER (ASSIGNED T w ) PART OF SURFACE 
TABLES OR UPPER (ASSIGNED B^) PART 
BY COMPARING CURRENT TEMPERATURE AT 
THIS STATION TO TEMPERATURES IN LOW- 
EST ASSIGNED B^ ENTRIES IN TABLES. 

IF CURRENT T w is GREATER, GO TO B' 
SECTION OF TABLES AND ITERATE ON §' 

IF CURRENT T w IS LESS, GO TO T w SEC- 
TION OF TABLES AND ITERATE ON T„. 

w 


ASSIGNED Bq SURFACE ENERGY BALANCE 
ITERATION 


t 

ASSIGNED T w 
ITERATION NOT 
SHOWN 




CONTINUED ON NEXT PAGE 


Figure 24. (continued) 
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BALANCED TO ACCEPTABLE TOLERANCE? 


i 


CORRECT B£, LIMITING CORRECTION TO TABLE 


INTERVAL WIDTHS, REPEAT TO CONVERGENCE 



ADJUST SIZE OF SURFACE NODAL BOX 

I 

ON LAST SURFACE POINT? 




INCREMENT COLUMN INDEX: I « I + 1 


RETURN TO MACABRE 


Figure 24 


(concluded) 




APPENDIX A 


TRANSPORT PROPERTIES FROM THE ACE PACKAGE 

When the ACE package of subroutines is used to generate transport 
property data for subsequent use in the transfer coefficient subroutine 
RUCH, the transport properties are calculated by ACE from expressions 
which derive from simple kinetic theory and the particular multicomponent 
diffusion representation previously discussed in Section 5.3.3. The 
development of these expressions is discussed in detail in Reference 19. 

A brief summary of this development, and the resulting expressions, are 
presented in this section. 

Diffusion Coefficients - In Section 5.3.3, a bifurcation approxima- 
tion for binary diffusion coefficients was mentioned which characterized 
multicomponent diffusion phenomena with reasonable accuracy without unduly 
complicating the system of equations to be solved. This simplification is 
achieved through a correlation for binary diffusion coefficients of the 
form 


M . . = 

il 


F . F v 
i 1 


(A-l) 


where D is a reference diffusion coefficient and the F^ are diffusion fac- 
tors. These quantities are discussed in detail in Reference 19. The 
incorporation of (A-l) in the S tef an-Maxwell relation for mass diffusion 
fluxes indicates that the diffusion flux of species i may be written in 
terms of only properties of species i and global system properties. Sub- 
ject to a few simplifying assumptions (Reference 19) , this expression 
for may be written 


pD)J. 
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where 


2. =», i x i /F i p 2 

(A- 3) 
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(A-4) 


w i - Z) x j F j 

j 

v 2 = S^ X j /F j (A-5) 

j 

The accuracy of this formulation is examined in Reference 19 for a variety of 
chemical systems. It is shown that the calculated by Equation (A-l) re- 

present a very substantial improvement over equal diffusion coefficients 
when compared to exact values calculat. 1 directly from kinetic theory. The 
calculation of the mixture viscosity and thermal conductivity is based 
on the diffusion factors given by (A-l) , and these will be discussed in the 
following paragraphs. 

Mixture Viscosity - The expression employed by the ACE program to 
calculate the mixture viscosity derives from rigorous first order kinetic 
theory (Reference 20) , subject to a few simplifying assumptions, as dis- 
cussed in Reference 19 . 


^mix 


I 


-E 

i=l 




RTy 

x i +1 ’ 385 




j=i 



(A-6) 


where is the viscosity of the pure species i. 


in terms of the self diffusion coefficients 


£ . . 
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The vk may be expressed 


= 


6A ii Pi,&ii 


(A-7) 


where At. is a ratio of collision integrals based on a Lennard-Jones inter- 
molecular potential. Substituting (A-l) and (A-7) into (A-6) results in the 
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following expression for the viscosity of the multicomponent mixture. 


mix 


= p£y 
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(A- 8) 


and this is the expression utilized to calculate the rixture viscosity out- 
put by the ACE program. 

Mixture Thermal Conductivity - The thermal conductivity in a poly- 
atomic gas mixture may be represented by (Reference 20) 


k . = k . + k. . (A-9) 

mix mono-mix int 

where k mono _ m -L X the thermal conductivity in a mixture computed neglecting 
all internal degrees of freedom and k^ nt is the contribution to the thermal 
conductivity of the mixture due to the internal degrees of freedom of the 
molecules. A simplified expression for the mono-mixture thermal conductivity 
can be derived in a manner similar to the procedure previously discussed for 
the mixture viscosity. This simplified expression is (from Reference 19) 
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mono-mix 
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x k 
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( A-10 ) 


where k. is the thermal conductivity of the pure species i neglecting 

all internal degrees of freedom of the molecule. The k^ mono ma Y be expressed 
in terms of the y^ as per 


i mono 


15 

T ^ 


“i 


(A— 11) 


The contribution to the thermal conductivity from the internal degrees of 
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freedom may be expressed as (from Reference 19) 


r In 5 r\ 

px i ~m ypi 2 

k int T 

1=1 V -L 

2L» (A-12) 

j=l 

By combining (A-l) with (A-9) through (A-12) , the mixture thermal conductivity 
may be written as 
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(A-13) 


where \x^ and y 2 are given by (A- 4) and (A-5) respectively, and 
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3 
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(A-14 ) 


( A- 1 5 ) 


Thus, (A-13) is the expression utilized to calculate the mixture thermal 
conductivity output by the ACE program. 

Also calculated and output by the ACE program are the Prandtl and 
Schmidt numbers which are defined here as 

Pr k C p-frozen 

Sc = — 

y 2 pD 

The transport properties calculated by the ACE program are all based 
on the bifurcation approximation for the binary diffusion coefficients ex- 
pressed in (A-l) . This is so even for closed system calculations (in which 
case diffusion phenomena need not be considered to calculate the chemical 
and thermodynamic state of the system) and for open system calculations for 


( A- 1 6 ) 
(A- 17) 
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which equal diffusion coefficients are assumed (Section 5-3.2). From the 
equations presented, it may be observed that the properties calculated are 
highly dependent on the diffusion factors, F^ . Three alternate methods for 
prescribing the F^ were discussed in Section 5.3.3. The use of the diffusion 
factor correlation (63) with resident values of and e (which were 

derived primarily from consideration of species diffusion coefficients) should 
result in reasonably accurate values of other transport properties . Alter- 
nately, the correlation (63) may be used with values of ^ re f and e derived 
by correlating available transport property data for the particular system 
of interest. If transport properties of maximum accuracy are to be cal- 
culated, then the diffusion factors should be input individually for each 
species. These data may be obtained from tabulations such as Reference 21. 
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